diff --git a/.mailmap b/.mailmap index 0ddd4f8ce..0adba676b 100644 --- a/.mailmap +++ b/.mailmap @@ -1,7 +1,6 @@ Maximilian Linhoff -Michele Peresano Michele Peresano -Michele Peresano Michele Peresano +Michele Peresano Michele Peresano Thomas Vuillaume vuillaut Thomas Vuillaume Thomas Vuillaume diff --git a/CHANGES.rst b/CHANGES.rst index 8294cebf4..90e01c2c7 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -1,3 +1,40 @@ +pyirf v0.11.0 (2024-05-14) +========================== + +Bug Fixes +--------- + +- Fix ``pyirf.benchmarks.energy_bias_resolution_from_energy_dispersion``. + This function was not adapted to the now correct normalization of the + energy dispersion matrix, resulting in wrong results on the now correct + matrices. [`#268 `__] + + +New Features +------------ + +- Adds an extrapolator for parametrized compontents utilizing blending over visible edges, resulting + in a smooth extrapolation compared to the NearestSimplexExtrapolator. [`#253 `__] + +- Ignore warnings about invalid floating point operations when calculating `n_signal` and `n_signal_weigthed` because the relative sensitivty is frequently NaN. [`#264 `__] + +- Add basic combinatoric fill-value handling for RAD_MAX estimation. [`#282 `__] + + +Maintenance +----------- + +- Clarified some documentation. [`#266 `__] + +- Add support for astropy 6.0. [`#271 `__] + +- Added filling of CREF keyword (IRF axis order) in the output files [`#275 `__] + + + +Refactoring and Optimization +---------------------------- + pyirf v0.10.1 (2023-09-15) ========================== diff --git a/README.rst b/README.rst index bdb35a58d..310e06a12 100644 --- a/README.rst +++ b/README.rst @@ -29,3 +29,20 @@ please cite it by using the corresponding DOI: - latest : |doilatest| - all versions: `Please visit Zenodo `_ and select the correct version + +At this point, our latest publication is the `2023 ICRC proceeding `_, which you can +cite using the following bibtex entry, especially if using functionalities from ``pyirf.interpolation``: + +.. code:: + + @inproceedings{pyirf-icrc-2023, + author = {Dominik, Rune Michael and Linhoff, Maximilian and Sitarek, Julian}, + title = {Interpolation of Instrument Response Functions for the Cherenkov Telescope Array in the Context of pyirf}, + usera = {for the CTA Consortium}, + doi = {10.22323/1.444.0703}, + booktitle = {Proceedings, 38th International Cosmic Ray Conference}, + year=2023, + volume={444}, + number={618}, + location={Nagoya, Japan}, + } diff --git a/docs/changes/253.feature.rst b/docs/changes/253.feature.rst deleted file mode 100644 index 6f10240a0..000000000 --- a/docs/changes/253.feature.rst +++ /dev/null @@ -1,2 +0,0 @@ -Adds an extrapolator for parametrized compontents utilizing blending over visible edges, resulting -in a smooth extrapolation compared to the NearestSimplexExtrapolator. diff --git a/docs/changes/264.feature.rst b/docs/changes/264.feature.rst deleted file mode 100644 index cde1ec34c..000000000 --- a/docs/changes/264.feature.rst +++ /dev/null @@ -1 +0,0 @@ -Ignore warnings about invalid floating point operations when calculating `n_signal` and `n_signal_weigthed` because the relative sensitivty is frequently NaN. diff --git a/docs/changes/266.maintenance.rst b/docs/changes/266.maintenance.rst deleted file mode 100644 index f5fdd08cc..000000000 --- a/docs/changes/266.maintenance.rst +++ /dev/null @@ -1 +0,0 @@ -Clarified some documentation. diff --git a/docs/changes/271.maintenance.rst b/docs/changes/271.maintenance.rst deleted file mode 100644 index 908e1b838..000000000 --- a/docs/changes/271.maintenance.rst +++ /dev/null @@ -1 +0,0 @@ -Add support for astropy 6.0. diff --git a/docs/changes/268.bugfix.rst b/docs/changes/273.bugfix.rst similarity index 63% rename from docs/changes/268.bugfix.rst rename to docs/changes/273.bugfix.rst index 7458041a1..607ffa71f 100644 --- a/docs/changes/268.bugfix.rst +++ b/docs/changes/273.bugfix.rst @@ -1,4 +1,4 @@ -Fix ``pyirf.benchmarks.energy_bias_resolution_from_energy_dispersion``. +Fix ``pyirf.irfs.energy_dispersion.energy_dispersion_to_migration``. This function was not adapted to the now correct normalization of the energy dispersion matrix, resulting in wrong results on the now correct -matrices. +matrices. \ No newline at end of file diff --git a/docs/changes/275.maintenance.rst b/docs/changes/275.maintenance.rst deleted file mode 100644 index 71ce09633..000000000 --- a/docs/changes/275.maintenance.rst +++ /dev/null @@ -1 +0,0 @@ -Added filling of CREF keyword (IRF axis order) in the output files diff --git a/docs/changes/276.feature.rst b/docs/changes/276.feature.rst new file mode 100644 index 000000000..83a61ea65 --- /dev/null +++ b/docs/changes/276.feature.rst @@ -0,0 +1 @@ +Add function to make 3d background for lon/lat coordinates. diff --git a/docs/changes/279.maintenance.rst b/docs/changes/279.maintenance.rst deleted file mode 100644 index b8a69c5d2..000000000 --- a/docs/changes/279.maintenance.rst +++ /dev/null @@ -1 +0,0 @@ -Temporarily fix scipy to <1.12 until gammapy supports API changes introduced with scipy 1.12. 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/docs/changes/289.maintenance.rst b/docs/changes/289.maintenance.rst new file mode 100644 index 000000000..379a7e43f --- /dev/null +++ b/docs/changes/289.maintenance.rst @@ -0,0 +1,3 @@ +Make pyirf compatible with numpy 2.0. + +Replace deprecated ``logging.warn`` with ``logging.warning``. diff --git a/docs/changes/290.api.rst b/docs/changes/290.api.rst new file mode 100644 index 000000000..882b2ba52 --- /dev/null +++ b/docs/changes/290.api.rst @@ -0,0 +1,3 @@ +Make it possible to pass multiple quantiles to ``pyirf.benchmarks.angular_resolution``, calculating all of them. + +The column name(s) in the output now include(s) the percentage value of the calculated quantile, e.g. ``angular_resolution_68``. diff --git a/docs/conf.py b/docs/conf.py index c8f4fc96d..59be578e7 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -84,8 +84,7 @@ html_theme = "sphinx_rtd_theme" html_theme_options = { - 'canonical_url': 'https://cta-observatory.github.io/pyirf', - 'display_version': True, + 'canonical_url': 'https://pyirf.readthedocs.io/', } intersphinx_mapping = { diff --git a/examples/comparison_with_EventDisplay.ipynb b/examples/comparison_with_EventDisplay.ipynb index 71c2d59c4..564b86379 100644 --- a/examples/comparison_with_EventDisplay.ipynb +++ b/examples/comparison_with_EventDisplay.ipynb @@ -528,7 +528,7 @@ "\n", "plt.errorbar(\n", " 0.5 * (ang_res['reco_energy_low'] + ang_res['reco_energy_high']).to_value(u.TeV),\n", - " ang_res['angular_resolution'].to_value(u.deg),\n", + " ang_res['angular_resolution_68'].to_value(u.deg),\n", " xerr=0.5 * (ang_res['reco_energy_high'] - ang_res['reco_energy_low']).to_value(u.TeV),\n", " ls='',\n", " label='pyirf'\n", diff --git a/pyirf/benchmarks/angular_resolution.py b/pyirf/benchmarks/angular_resolution.py index 2f3445cfe..d87f60ba9 100644 --- a/pyirf/benchmarks/angular_resolution.py +++ b/pyirf/benchmarks/angular_resolution.py @@ -1,3 +1,4 @@ +from collections.abc import Sequence import numpy as np from astropy.table import QTable from scipy.stats import norm @@ -10,7 +11,8 @@ def angular_resolution( - events, energy_bins, + events, + energy_bins, energy_type="true", quantile=ONE_SIGMA_QUANTILE, ): @@ -20,6 +22,8 @@ def angular_resolution( This implementation corresponds to the 68% containment of the angular distance distribution. + Passing a list of quantiles results in all the quantiles being calculated. + Parameters ---------- events : astropy.table.QTable @@ -29,8 +33,8 @@ def angular_resolution( energy_type: str Either "true" or "reco" energy. Default is "true". - quantile : float - Which quantile to use for the angular resolution, + quantile : list(float) + Which quantile(s) to use for the angular resolution, by default, the containment of the 1-sigma region of the normal distribution (~68%) is used. @@ -52,8 +56,13 @@ def angular_resolution( result[f"{energy_key}_center"] = 0.5 * (energy_bins[:-1] + energy_bins[1:]) result["n_events"] = 0 - key = "angular_resolution" - result[key] = u.Quantity(np.nan, table["theta"].unit) + if not isinstance(quantile, Sequence): + quantile = [quantile] + + keys = [f"angular_resolution_{value * 100:.0f}" for value in quantile] + + for key in keys: + result[key] = u.Quantity(np.nan, table["theta"].unit) # if we get an empty input (no selected events available) # we return the table filled with NaNs @@ -63,6 +72,9 @@ def angular_resolution( # use groupby operations to calculate the percentile in each bin by_bin = table[valid].group_by(bin_index[valid]) for bin_idx, group in zip(by_bin.groups.keys, by_bin.groups): - result[key][bin_idx] = np.nanquantile(group["theta"], quantile) result["n_events"][bin_idx] = len(group) + quantile_values = np.nanquantile(group["theta"], quantile) + for key, value in zip(keys, quantile_values): + result[key][bin_idx] = value + return result diff --git a/pyirf/benchmarks/energy_bias_resolution.py b/pyirf/benchmarks/energy_bias_resolution.py index 675c9af18..4479549c3 100644 --- a/pyirf/benchmarks/energy_bias_resolution.py +++ b/pyirf/benchmarks/energy_bias_resolution.py @@ -4,6 +4,7 @@ import astropy.units as u from ..binning import calculate_bin_indices, UNDERFLOW_INDEX, OVERFLOW_INDEX +from ..compat import COPY_IF_NEEDED NORM_LOWER_SIGMA, NORM_UPPER_SIGMA = norm(0, 1).cdf([-1, 1]) @@ -83,8 +84,10 @@ def energy_bias_resolution( """ # create a table to make use of groupby operations - table = QTable(events[["true_energy", "reco_energy"]], copy=False) - table["rel_error"] = (events["reco_energy"] / events["true_energy"]).to_value(u.one) - 1 + table = QTable(events[["true_energy", "reco_energy"]], copy=COPY_IF_NEEDED) + table["rel_error"] = (events["reco_energy"] / events["true_energy"]).to_value( + u.one + ) - 1 energy_key = f"{energy_type}_energy" @@ -102,7 +105,6 @@ def energy_bias_resolution( # we return the table filled with NaNs return result - # use groupby operations to calculate the percentile in each bin bin_index, valid = calculate_bin_indices(table[energy_key], energy_bins) by_bin = table.group_by(bin_index) @@ -115,6 +117,7 @@ def energy_bias_resolution( result["resolution"][bin_idx] = resolution_function(group["rel_error"]) return result + def energy_bias_resolution_from_energy_dispersion( energy_dispersion, migration_bins, @@ -147,7 +150,7 @@ def energy_bias_resolution_from_energy_dispersion( low, median, high = np.interp( [NORM_LOWER_SIGMA, MEDIAN, NORM_UPPER_SIGMA], cdf[energy_bin, :, fov_bin], - migration_bins[1:] # cdf is defined at upper bin edge + migration_bins[1:], # cdf is defined at upper bin edge ) bias[energy_bin, fov_bin] = median - 1 resolution[energy_bin, fov_bin] = 0.5 * (high - low) diff --git a/pyirf/benchmarks/tests/test_angular_resolution.py b/pyirf/benchmarks/tests/test_angular_resolution.py index 269c6c46b..e66ed50dd 100644 --- a/pyirf/benchmarks/tests/test_angular_resolution.py +++ b/pyirf/benchmarks/tests/test_angular_resolution.py @@ -8,17 +8,17 @@ def test_empty_angular_resolution(): from pyirf.benchmarks import angular_resolution - events = QTable({ - 'true_energy': [] * u.TeV, - 'theta': [] * u.deg, - }) - - table = angular_resolution( - events, - [1, 10, 100] * u.TeV + events = QTable( + { + "true_energy": [] * u.TeV, + "theta": [] * u.deg, + } ) - assert np.all(np.isnan(table["angular_resolution"])) + table = angular_resolution(events, [1, 10, 100] * u.TeV) + + assert np.all(np.isnan(table["angular_resolution_68"])) + @pytest.mark.parametrize("unit", (u.deg, u.rad)) def test_angular_resolution(unit): @@ -32,20 +32,24 @@ def test_angular_resolution(unit): true_resolution = np.append(np.full(N, TRUE_RES_1), np.full(N, TRUE_RES_2)) - rng = np.random.default_rng(0) - events = QTable({ - 'true_energy': np.concatenate([ - [0.5], # below bin 1 to test with underflow - np.full(N - 1, 5.0), - np.full(N - 1, 50.0), - [500], # above bin 2 to test with overflow - ]) * u.TeV, - 'theta': np.abs(rng.normal(0, true_resolution)) * u.deg - }) + events = QTable( + { + "true_energy": np.concatenate( + [ + [0.5], # below bin 1 to test with underflow + np.full(N - 1, 5.0), + np.full(N - 1, 50.0), + [500], # above bin 2 to test with overflow + ] + ) + * u.TeV, + "theta": np.abs(rng.normal(0, true_resolution)) * u.deg, + } + ) - events['theta'] = events['theta'].to(unit) + events["theta"] = events["theta"].to(unit) # add nans to test if nans are ignored events["true_energy"].value[N // 2] = np.nan @@ -53,7 +57,7 @@ def test_angular_resolution(unit): bins = [1, 10, 100] * u.TeV table = angular_resolution(events, bins) - ang_res = table['angular_resolution'].to(u.deg) + ang_res = table["angular_resolution_68"].to(u.deg) assert len(ang_res) == 2 assert u.isclose(ang_res[0], TRUE_RES_1 * u.deg, rtol=0.05) assert u.isclose(ang_res[1], TRUE_RES_2 * u.deg, rtol=0.05) @@ -64,9 +68,26 @@ def test_angular_resolution(unit): # 2 sigma coverage interval quantile = norm(0, 1).cdf(2) - norm(0, 1).cdf(-2) table = angular_resolution(events, bins, quantile=quantile) - ang_res = table['angular_resolution'].to(u.deg) - + ang_res = table["angular_resolution_95"].to(u.deg) assert len(ang_res) == 2 - assert u.isclose(ang_res[0], 2 * TRUE_RES_1 * u.deg, rtol=0.05) assert u.isclose(ang_res[1], 2 * TRUE_RES_2 * u.deg, rtol=0.05) + + # 25%, 50%, 90% coverage interval + table = angular_resolution(events, bins, quantile=[0.25, 0.5, 0.9]) + cov_25 = table["angular_resolution_25"].to(u.deg) + assert len(cov_25) == 2 + assert u.isclose( + cov_25[0], norm(0, TRUE_RES_1).interval(0.25)[1] * u.deg, rtol=0.05 + ) + assert u.isclose( + cov_25[1], norm(0, TRUE_RES_2).interval(0.25)[1] * u.deg, rtol=0.05 + ) + + cov_50 = table["angular_resolution_50"].to(u.deg) + assert u.isclose(cov_50[0], norm(0, TRUE_RES_1).interval(0.5)[1] * u.deg, rtol=0.05) + assert u.isclose(cov_50[1], norm(0, TRUE_RES_2).interval(0.5)[1] * u.deg, rtol=0.05) + + cov_90 = table["angular_resolution_90"].to(u.deg) + assert u.isclose(cov_90[0], norm(0, TRUE_RES_1).interval(0.9)[1] * u.deg, rtol=0.05) + assert u.isclose(cov_90[1], norm(0, TRUE_RES_2).interval(0.9)[1] * u.deg, rtol=0.05) diff --git a/pyirf/binning.py b/pyirf/binning.py index 1aa3ecd02..0eb4413c4 100644 --- a/pyirf/binning.py +++ b/pyirf/binning.py @@ -7,6 +7,8 @@ import astropy.units as u from astropy.table import QTable +from .compat import COPY_IF_NEEDED + #: Index returned by `calculate_bin_indices` for underflown values UNDERFLOW_INDEX = np.iinfo(np.int64).min @@ -37,12 +39,12 @@ def join_bin_lo_hi(bin_lo, bin_hi): The joint bins """ - if np.allclose(bin_lo[...,1:], bin_hi[...,:-1], rtol=1.e-5): - last_axis=len(bin_lo.shape)-1 - bins = np.concatenate((bin_lo, bin_hi[...,-1:]), axis=last_axis) + if np.allclose(bin_lo[..., 1:], bin_hi[..., :-1], rtol=1.0e-5): + last_axis = len(bin_lo.shape) - 1 + bins = np.concatenate((bin_lo, bin_hi[..., -1:]), axis=last_axis) return bins else: - raise ValueError('Not matching bin edges') + raise ValueError("Not matching bin edges") def split_bin_lo_hi(bins): @@ -62,10 +64,11 @@ def split_bin_lo_hi(bins): bin_hi: np.array or u.Quantity Hi bin edges array """ - bin_lo=bins[...,:-1] - bin_hi=bins[...,1:] + bin_lo = bins[..., :-1] + bin_hi = bins[..., 1:] return bin_lo, bin_hi + def add_overflow_bins(bins, positive=True): """ Add under and overflow bins to a bin array. @@ -124,7 +127,7 @@ def create_bins_per_decade(e_min, e_max, bins_per_decade=5): # include endpoint if reasonably close eps = step / 10000 bins = 10 ** np.arange(log_lower, log_upper + eps, step) - return u.Quantity(bins, e_min.unit, copy=False) + return u.Quantity(bins, e_min.unit, copy=COPY_IF_NEEDED) def calculate_bin_indices(data, bins): @@ -156,7 +159,7 @@ def calculate_bin_indices(data, bins): valid: np.ndarray[bool] Boolean mask indicating if a given value belongs into one of the defined bins. - False indicates that an entry fell into the over- or underflow bins. + False indicates that an entry fell into the over- or underflow bins. """ if hasattr(data, "unit"): @@ -169,8 +172,8 @@ def calculate_bin_indices(data, bins): n_bins = len(bins) - 1 idx = np.digitize(data, bins) - 1 - underflow = (idx < 0) - overflow = (idx >= n_bins) + underflow = idx < 0 + overflow = idx >= n_bins idx[underflow] = UNDERFLOW_INDEX idx[overflow] = OVERFLOW_INDEX valid = ~underflow & ~overflow @@ -211,11 +214,11 @@ def create_histogram_table(events, bins, key="reco_energy"): hist["n_weighted"] = hist["n"] # create counts per particle type, only works if there is at least 1 event - if 'particle_type' in events.colnames and len(events) > 0: - by_particle = events.group_by('particle_type') + if "particle_type" in events.colnames and len(events) > 0: + by_particle = events.group_by("particle_type") for group_key, group in zip(by_particle.groups.keys, by_particle.groups): - particle = group_key['particle_type'] + particle = group_key["particle_type"] hist["n_" + particle], _ = np.histogram(group[key], bins) @@ -282,8 +285,11 @@ def resample_histogram1d(data, old_edges, new_edges, axis=0): cumsum /= norm f_integral = interp1d( - old_edges, cumsum, bounds_error=False, - fill_value=(0, 1), kind="quadratic", + old_edges, + cumsum, + bounds_error=False, + fill_value=(0, 1), + kind="quadratic", axis=axis, ) diff --git a/pyirf/compat.py b/pyirf/compat.py new file mode 100644 index 000000000..a52cd60a5 --- /dev/null +++ b/pyirf/compat.py @@ -0,0 +1,10 @@ +import numpy as np +from packaging.version import Version + + +# in numpy 1.x, copy=False allows copying if it cannot be avoided +# in numpy 2.0, copy=False raises an error when the copy cannot be avoided +# copy=None is a new option in numpy 2.0 for the previous behavior of copy=False +COPY_IF_NEEDED = None +if Version(np.__version__) < Version("2.0.0.dev"): + COPY_IF_NEEDED = False 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/cuts.py b/pyirf/cuts.py index 4f83268b1..ffda0e906 100644 --- a/pyirf/cuts.py +++ b/pyirf/cuts.py @@ -4,17 +4,25 @@ import astropy.units as u from .binning import calculate_bin_indices, bin_center +from .compat import COPY_IF_NEEDED __all__ = [ - 'calculate_percentile_cut', - 'evaluate_binned_cut', - 'compare_irf_cuts', + "calculate_percentile_cut", + "evaluate_binned_cut", + "compare_irf_cuts", ] def calculate_percentile_cut( - values, bin_values, bins, fill_value, percentile=68, min_value=None, max_value=None, - smoothing=None, min_events=10, + values, + bin_values, + bins, + fill_value, + percentile=68, + min_value=None, + max_value=None, + smoothing=None, + min_events=10, ): """ Calculate cuts as the percentile of a given quantity in bins of another @@ -46,7 +54,7 @@ def calculate_percentile_cut( """ # create a table to make use of groupby operations # we use a normal table here to avoid astropy/astropy#13840 - table = Table({"values": values}, copy=False) + table = Table({"values": values}, copy=COPY_IF_NEEDED) unit = table["values"].unit # make sure units match @@ -84,10 +92,10 @@ def calculate_percentile_cut( cut_table["cut"].value[bin_idx] = value if smoothing is not None: - cut_table['cut'].value[:] = gaussian_filter1d( + cut_table["cut"].value[:] = gaussian_filter1d( cut_table["cut"].value, smoothing, - mode='nearest', + mode="nearest", ) return cut_table @@ -126,7 +134,7 @@ def evaluate_binned_cut(values, bin_values, cut_table, op): passes the bin specific cut given in cut table. """ if not isinstance(cut_table, QTable): - raise ValueError('cut_table needs to be an astropy.table.QTable') + raise ValueError("cut_table needs to be an astropy.table.QTable") bins = np.append(cut_table["low"], cut_table["high"][-1]) bin_index, valid = calculate_bin_indices(bin_values, bins) diff --git a/pyirf/interpolation/component_estimators.py b/pyirf/interpolation/component_estimators.py index 740180e60..06e123f37 100644 --- a/pyirf/interpolation/component_estimators.py +++ b/pyirf/interpolation/component_estimators.py @@ -1,19 +1,23 @@ """Classes to estimate (interpolate/extrapolate) actual IRF HDUs""" + import warnings import astropy.units as u import numpy as np -from pyirf.utils import cone_solid_angle from scipy.spatial import Delaunay + from .base_extrapolators import DiscretePDFExtrapolator, ParametrizedExtrapolator from .base_interpolators import ( DiscretePDFInterpolator, ParametrizedInterpolator, PDFNormalization, ) +from ..compat import COPY_IF_NEEDED from .griddata_interpolator import GridDataInterpolator +from .nearest_neighbor_searcher import BaseNearestNeighborSearcher from .quantile_interpolator import QuantileInterpolator +from .utils import find_nearest_simplex __all__ = [ "BaseComponentEstimator", @@ -427,9 +431,9 @@ def __init__( self.min_effective_area = min_effective_area # remove zeros and log it - effective_area[ - effective_area < self.min_effective_area - ] = self.min_effective_area + effective_area[effective_area < self.min_effective_area] = ( + self.min_effective_area + ) effective_area = np.log(effective_area) super().__init__( @@ -467,7 +471,7 @@ def __call__(self, target_point): # remove entries manipulated by min_effective_area aeff_interp[aeff_interp < self.min_effective_area] = 0 - return u.Quantity(aeff_interp, u.m**2, copy=False) + return u.Quantity(aeff_interp, u.m**2, copy=COPY_IF_NEEDED) class RadMaxEstimator(ParametrizedComponentEstimator): @@ -477,6 +481,7 @@ def __init__( self, grid_points, rad_max, + fill_value=None, interpolator_cls=GridDataInterpolator, interpolator_kwargs=None, extrapolator_cls=None, @@ -495,6 +500,13 @@ def __init__( Grid of theta cuts. Dimensions but the first can in principle be freely chosen. Class is RAD_MAX_2D compatible, which would require shape=(n_points, n_energy_bins, n_fov_offset_bins). + fill_val: + Indicator of fill-value handling. If None, fill values are regarded as + normal values and no special handeling is applied. If "infer", fill values + will be infered as max of all values, if a value is provided, + it is used to flag fill values. Flagged fill-values are + not used for interpolation. Fill-value handling is only supported in + up to two grid dimensions. interpolator_cls: pyirf interpolator class, defaults to GridDataInterpolator. interpolator_kwargs: dict @@ -523,6 +535,26 @@ def __init__( extrapolator_kwargs=extrapolator_kwargs, ) + self.params = rad_max + + if fill_value is None: + self.fill_val = None + elif fill_value == "infer": + self.fill_val = np.max(self.params) + else: + self.fill_val = fill_value + + # Raise error if fill-values should be handled in >=3 dims + if self.fill_val and self.grid_dim >= 3: + raise ValueError( + "Fill-value handling only supported in up to two grid dimensions." + ) + + # If fill-values should be handled in 2D, construct a trinangulation + # to later determine in which simplex the target values is + if self.fill_val and (self.grid_dim == 2): + self.triangulation = Delaunay(self.grid_points) + def __call__(self, target_point): """ Estimating rad max table at target_point, inter-/extrapolates as needed and @@ -540,8 +572,99 @@ def __call__(self, target_point): effective areas. For RAD_MAX_2D of shape (n_energy_bins, n_fov_offset_bins) """ + # First, construct estimation without handling fill-values + full_estimation = super().__call__(target_point) + # Safeguard against extreme extrapolation cases + np.clip(full_estimation, 0, None, out=full_estimation) + + # Early exit if fill_values should not be handled + if not self.fill_val: + return full_estimation + + # Early exit if a nearest neighbor estimation would be overwritten + # Complex setup is needed to catch settings where the user mixes approaches and + # e.g. uses nearest neighbors for extrapolation and an actual interpolation otherwise + if self.grid_dim == 1: + if ( + (target_point < self.grid_points.min()) + or (target_point > self.grid_points.max()) + ) and issubclass(self.extrapolator.__class__, BaseNearestNeighborSearcher): + return full_estimation + elif issubclass(self.interpolator.__class__, BaseNearestNeighborSearcher): + return full_estimation + elif self.grid_dim == 2: + target_simplex = self.triangulation.find_simplex(target_point) + + if (target_simplex == -1) and issubclass( + self.extrapolator.__class__, BaseNearestNeighborSearcher + ): + return full_estimation + elif issubclass(self.interpolator.__class__, BaseNearestNeighborSearcher): + return full_estimation + + # Actual fill-value handling + if self.grid_dim == 1: + # Locate target in grid + if target_point < self.grid_points.min(): + segment_inds = np.array([0, 1], "int") + elif target_point > self.grid_points.max(): + segment_inds = np.array([-2, -1], "int") + else: + target_bin = np.digitize( + target_point.squeeze(), self.grid_points.squeeze() + ) + segment_inds = np.array([target_bin - 1, target_bin], "int") + + mask_left = self.params[segment_inds[0]] >= self.fill_val + mask_right = self.params[segment_inds[1]] >= self.fill_val + # Indicate, wether one of the neighboring entries is a fill-value + mask = np.logical_or(mask_left, mask_right) + elif self.grid_dim == 2: + # Locate target + target_simplex = self.triangulation.find_simplex(target_point) + + if target_simplex == -1: + target_simplex = find_nearest_simplex(self.triangulation, target_point) + + simplex_nodes_indices = self.triangulation.simplices[ + target_simplex + ].squeeze() + + mask0 = self.params[simplex_nodes_indices[0]] >= self.fill_val + mask1 = self.params[simplex_nodes_indices[1]] >= self.fill_val + mask2 = self.params[simplex_nodes_indices[2]] >= self.fill_val + + # This collected mask now counts for each entry in the estimation how many + # of the entries used for extrapolation contained fill-values + intermediate_mask = ( + mask0.astype("int") + mask1.astype("int") + mask2.astype("int") + ) + mask = np.full_like(intermediate_mask, True, dtype=bool) + + # Simplest cases: All or none entries were fill-values, so either return + # a fill-value or the actual estimation + mask[intermediate_mask == 0] = False + mask[intermediate_mask == 3] = True + + # If two out of three values were fill-values return a fill-value as estimate + mask[intermediate_mask == 2] = True + + # If only one out of three values was a fill-value use the smallest value of the + # remaining two + mask[intermediate_mask == 1] = False + full_estimation = np.where( + intermediate_mask[np.newaxis, :] == 1, + np.min(self.params[simplex_nodes_indices], axis=0), + full_estimation, + ) + + # Set all flagged values to fill-value + full_estimation[mask[np.newaxis, :]] = self.fill_val + + # Safeguard against extreme extrapolation cases + full_estimation[full_estimation > self.fill_val] = self.fill_val - return super().__call__(target_point) + return full_estimation class EnergyDispersionEstimator(DiscretePDFComponentEstimator): @@ -732,5 +855,7 @@ def __call__(self, target_point): # Undo normalisation to get a proper PSF and return return u.Quantity( - np.swapaxes(interpolated_psf_normed, -1, self.axis), u.sr**-1, copy=False + np.swapaxes(interpolated_psf_normed, -1, self.axis), + u.sr**-1, + copy=COPY_IF_NEEDED, ) diff --git a/pyirf/interpolation/tests/test_component_estimator_specific_classes.py b/pyirf/interpolation/tests/test_component_estimator_specific_classes.py index c1f77913b..af786c4e6 100644 --- a/pyirf/interpolation/tests/test_component_estimator_specific_classes.py +++ b/pyirf/interpolation/tests/test_component_estimator_specific_classes.py @@ -1,8 +1,10 @@ import astropy.units as u import numpy as np -from pyirf.utils import cone_solid_angle +import pytest from scipy.stats import expon +from pyirf.utils import cone_solid_angle + def test_EnergyDispersionEstimator(prod5_irfs): from pyirf.interpolation import EnergyDispersionEstimator, QuantileInterpolator @@ -74,7 +76,9 @@ def hist(pnt): interp = estimator(target_point=zen_pnt[[1]]) - probability = (interp * omegas[np.newaxis, np.newaxis, np.newaxis, ...]).to_value(u.one) + probability = (interp * omegas[np.newaxis, np.newaxis, np.newaxis, ...]).to_value( + u.one + ) assert np.max(probability) <= 1 assert np.min(probability) >= 0 @@ -177,6 +181,7 @@ def test_RadMaxEstimator(): estimator = RadMaxEstimator( grid_points=grid_points, rad_max=rad_max, + fill_value=None, interpolator_cls=GridDataInterpolator, interpolator_kwargs=None, extrapolator_cls=None, @@ -186,3 +191,183 @@ def test_RadMaxEstimator(): assert interp.shape == (1, *rad_max_1.shape) assert np.allclose(interp, 1.5 * rad_max_1) + + +def test_RadMaxEstimator_fill_val_handling_1D(): + from pyirf.interpolation import ( + GridDataInterpolator, + ParametrizedNearestNeighborSearcher, + ParametrizedNearestSimplexExtrapolator, + RadMaxEstimator, + ) + + grid_points_1D = np.array([[0], [1], [2]]) + + rad_max_1 = np.array([[0.95, 0.95, 0.5, 0.95, 0.95], [0.95, 0.5, 0.3, 0.5, 0.95]]) + rad_max_2 = np.array([[0.95, 0.5, 0.3, 0.5, 0.95], [0.5, 0.3, 0.2, 0.9, 0.5]]) + rad_max_3 = np.array([[0.95, 0.4, 0.2, 0.4, 0.5], [0.5, 0.3, 0, 0.94, 0.6]]) + + rad_max_1D = np.array([rad_max_1, rad_max_2, rad_max_3]) + + truth_0 = np.array([[0.95, 0.95, 0.7, 0.95, 0.95], [0.95, 0.7, 0.4, 0.1, 0.95]]) + truth_1_5 = np.array([[0.95, 0.95, 0.4, 0.95, 0.95], [0.95, 0.4, 0.25, 0.7, 0.95]]) + + truth_4 = np.array([[0.95, 0.3, 0.1, 0.3, 0.95], [0.5, 0.3, 0, 0.95, 0.7]]) + + # State fill value + estim = RadMaxEstimator( + grid_points=grid_points_1D, + rad_max=rad_max_1D, + fill_value=0.95, + interpolator_cls=GridDataInterpolator, + interpolator_kwargs=None, + extrapolator_cls=ParametrizedNearestSimplexExtrapolator, + extrapolator_kwargs=None, + ) + + assert np.allclose(estim(np.array([-1])), truth_0) + assert np.allclose(estim(np.array([0.5])), truth_1_5) + assert np.allclose(estim(np.array([3])), truth_4) + + # Infer fill-val as max of rad-max vals + estim = RadMaxEstimator( + grid_points=grid_points_1D, + rad_max=rad_max_1D, + fill_value="infer", + interpolator_cls=GridDataInterpolator, + interpolator_kwargs=None, + extrapolator_cls=ParametrizedNearestSimplexExtrapolator, + extrapolator_kwargs=None, + ) + + assert np.allclose(estim(np.array([0.5])), truth_1_5) + assert np.allclose(estim(np.array([3])), truth_4) + + # Nearest neighbor cases + estim = RadMaxEstimator( + grid_points=grid_points_1D, + rad_max=rad_max_1D, + fill_value="infer", + interpolator_cls=ParametrizedNearestNeighborSearcher, + interpolator_kwargs=None, + extrapolator_cls=ParametrizedNearestNeighborSearcher, + extrapolator_kwargs=None, + ) + + assert np.allclose(estim(np.array([0.25])), rad_max_1) + assert np.allclose(estim(np.array([3])), rad_max_3) + + # Ignore fill values + estim = RadMaxEstimator( + grid_points=grid_points_1D, + rad_max=rad_max_1D, + fill_value=None, + interpolator_cls=GridDataInterpolator, + interpolator_kwargs=None, + extrapolator_cls=None, + extrapolator_kwargs=None, + ) + + assert np.allclose(estim(np.array([0.5])), (rad_max_1 + rad_max_2) / 2) + + +def test_RadMaxEstimator_fill_val_handling_2D(): + from pyirf.interpolation import ( + GridDataInterpolator, + ParametrizedNearestNeighborSearcher, + ParametrizedNearestSimplexExtrapolator, + RadMaxEstimator, + ) + + grid_points_2D = np.array([[0, 0], [1, 0], [0, 1]]) + + rad_max_1 = np.array([[0.95, 0.95, 0.5, 0.95, 0.95], [0.5, 0.5, 0.3, 0.5, 0.5]]) + rad_max_2 = np.array([[0.95, 0.95, 0.5, 0.5, 0.95], [0.95, 0.95, 0.95, 0.5, 0.95]]) + rad_max_3 = np.array([[0.95, 0.5, 0.5, 0.4, 0.5], [0.4, 0.95, 0, 0.5, 0.95]]) + + rad_max_2D = np.array([rad_max_1, rad_max_2, rad_max_3]) + + # Only test for combinatoric cases, thus inter- and extrapolation have the same + # result in this special test case. Correct estimation is checked elsewhere + truth = np.array([[0.95, 0.95, 0.5, 0.4, 0.95], [0.4, 0.95, 0, 0.5, 0.95]]) + + # State fill-value + estim = RadMaxEstimator( + grid_points=grid_points_2D, + rad_max=rad_max_2D, + fill_value=0.95, + interpolator_cls=GridDataInterpolator, + interpolator_kwargs=None, + extrapolator_cls=ParametrizedNearestSimplexExtrapolator, + extrapolator_kwargs=None, + ) + + assert np.allclose(estim(np.array([0.5, 0.5])), truth) + assert np.allclose(estim(np.array([-1, -1])), truth) + + # Infer fill-val as max of rad-max vals + estim = RadMaxEstimator( + grid_points=grid_points_2D, + rad_max=rad_max_2D, + fill_value="infer", + interpolator_cls=GridDataInterpolator, + interpolator_kwargs=None, + extrapolator_cls=ParametrizedNearestSimplexExtrapolator, + extrapolator_kwargs=None, + ) + + assert np.allclose(estim(np.array([0.5, 0.5])), truth) + assert np.allclose(estim(np.array([-1, -1])), truth) + + # Nearest neighbor cases + estim = RadMaxEstimator( + grid_points=grid_points_2D, + rad_max=rad_max_2D, + fill_value="infer", + interpolator_cls=ParametrizedNearestNeighborSearcher, + interpolator_kwargs=None, + extrapolator_cls=ParametrizedNearestNeighborSearcher, + extrapolator_kwargs=None, + ) + + assert np.allclose(estim(np.array([0.25, 0.25])), rad_max_1) + assert np.allclose(estim(np.array([0, 1.1])), rad_max_3) + + # Ignore fill-values + estim = RadMaxEstimator( + grid_points=grid_points_2D, + rad_max=rad_max_2D, + fill_value=None, + interpolator_cls=GridDataInterpolator, + interpolator_kwargs=None, + extrapolator_cls=None, + extrapolator_kwargs=None, + ) + + truth_interpolator = GridDataInterpolator(grid_points_2D, rad_max_2D) + + assert np.allclose( + estim(np.array([0.25, 0.25])), truth_interpolator(np.array([0.25, 0.25])) + ) + + +def test_RadMaxEstimator_fill_val_handling_3D(): + from pyirf.interpolation import GridDataInterpolator, RadMaxEstimator + + grid_points_3D = np.array([[0, 0, 0], [1, 0, 0], [0, 1, 0], [0, 0, 1]]) + + rad_max = np.array([[0.95], [0.95], [0.95], [0.95]]) + + with pytest.raises( + ValueError, + match="Fill-value handling only supported in up to two grid dimensions.", + ): + RadMaxEstimator( + grid_points=grid_points_3D, + rad_max=rad_max, + fill_value=0.95, + interpolator_cls=GridDataInterpolator, + interpolator_kwargs=None, + extrapolator_cls=None, + extrapolator_kwargs=None, + ) diff --git a/pyirf/io/__init__.py b/pyirf/io/__init__.py index aeefad535..a1e67c6f4 100644 --- a/pyirf/io/__init__.py +++ b/pyirf/io/__init__.py @@ -5,6 +5,7 @@ create_psf_table_hdu, create_rad_max_hdu, create_background_2d_hdu, + create_background_3d_hdu, ) @@ -16,4 +17,5 @@ "create_psf_table_hdu", "create_rad_max_hdu", "create_background_2d_hdu", + "create_background_3d_hdu", ] diff --git a/pyirf/io/gadf.py b/pyirf/io/gadf.py index e14c20670..764fc8118 100644 --- a/pyirf/io/gadf.py +++ b/pyirf/io/gadf.py @@ -260,6 +260,118 @@ def create_background_2d_hdu( return BinTableHDU(bkg, header=header, name=extname) +@u.quantity_input( + background=GADF_BACKGROUND_UNIT, reco_energy_bins=u.TeV, fov_offset_bins=u.deg, +) +def create_background_3d_hdu( + background_3d, + reco_energy_bins, + fov_lon_bins, + fov_lat_bins, + alignment="ALTAZ", + extname="BACKGROUND", + **header_cards, +): + """ + Create a fits binary table HDU in GADF format for the background 3d table. Assumes ALTAZ coordinates. + See the specification at + https://gamma-astro-data-formats.readthedocs.io/en/latest/irfs/full_enclosure/bkg/index.html#bkg-2d + + Parameters + ---------- + background_3d: astropy.units.Quantity[(MeV s sr)^-1] + Background rate, must have shape + (n_energy_bins, n_fov_offset_bins, n_fov_offset_bins) + reco_energy_bins: astropy.units.Quantity[energy] + Bin edges in reconstructed energy + fov_lon_bins: astropy.units.Quantity[angle] + Bin edges in the field of view system, becomes the DETX values + fov_lat_bins: astropy.units.Quantity[angle] + Bin edges in the field of view system, becomes the DETY values + alignment: str + Whether the FOV coordinates are aligned with the ALTAZ or RADEC system, more details at + https://gamma-astro-data-formats.readthedocs.io/en/latest/general/coordinates.html + extname: str + Name for BinTableHDU + **header_cards + Additional metadata to add to the header, use this to set e.g. TELESCOP or + INSTRUME. + """ + + bkg = QTable() + bkg["ENERG_LO"], bkg["ENERG_HI"] = binning.split_bin_lo_hi(reco_energy_bins[np.newaxis, :].to(u.TeV)) + bkg["DETX_LO"], bkg["DETX_HI"] = binning.split_bin_lo_hi(fov_lon_bins[np.newaxis, :].to(u.deg)) + bkg["DETY_LO"], bkg["DETY_HI"] = binning.split_bin_lo_hi(fov_lat_bins[np.newaxis, :].to(u.deg)) + # transpose as FITS uses opposite dimension order + bkg["BKG"] = background_3d.T[np.newaxis, ...].to(GADF_BACKGROUND_UNIT) + + # required header keywords + header = DEFAULT_HEADER.copy() + header["HDUCLAS1"] = "RESPONSE" + header["HDUCLAS2"] = "BKG" + header["HDUCLAS3"] = "FULL-ENCLOSURE" + header["HDUCLAS4"] = "BKG_3D" + if alignment in ["ALTAZ", "RADEC"]: + header["FOVALIGN"] = alignment + else: + raise ValueError(f"'alignment' must be one of 'ALTAZ' or 'RADEC', got {alignment}") + header["DATE"] = Time.now().utc.iso + _add_header_cards(header, **header_cards) + + return BinTableHDU(bkg, header=header, name=extname) + + +@u.quantity_input( + background=GADF_BACKGROUND_UNIT, reco_energy_bins=u.TeV, fov_offset_bins=u.deg, +) +def create_background_3d_hdu( + background_3d, + reco_energy_bins, + fov_offset_bins, + extname="BACKGROUND", + **header_cards, +): + """ + Create a fits binary table HDU in GADF format for the background 3d table. Assumes ALTAZ coordinates. + See the specification at + https://gamma-astro-data-formats.readthedocs.io/en/latest/irfs/full_enclosure/bkg/index.html#bkg-2d + + Parameters + ---------- + background_3d: astropy.units.Quantity[(MeV s sr)^-1] + Background rate, must have shape + (n_energy_bins, n_fov_offset_bins, n_fov_offset_bins) + reco_energy_bins: astropy.units.Quantity[energy] + Bin edges in reconstructed energy + fov_offset_bins: astropy.units.Quantity[angle] + Bin edges in the field of view offset. + extname: str + Name for BinTableHDU + **header_cards + Additional metadata to add to the header, use this to set e.g. TELESCOP or + INSTRUME. + """ + + bkg = QTable() + bkg["ENERG_LO"], bkg["ENERG_HI"] = binning.split_bin_lo_hi(reco_energy_bins[np.newaxis, :].to(u.TeV)) + bkg["DETX_LO"], bkg["DETX_HI"] = binning.split_bin_lo_hi(fov_offset_bins[np.newaxis, :].to(u.deg)) + bkg["DETY_LO"], bkg["DETY_HI"] = binning.split_bin_lo_hi(fov_offset_bins[np.newaxis, :].to(u.deg)) + # transpose as FITS uses opposite dimension order + bkg["BKG"] = background_3d.T[np.newaxis, ...].to(GADF_BACKGROUND_UNIT) + + # required header keywords + header = DEFAULT_HEADER.copy() + header["HDUCLAS1"] = "RESPONSE" + header["HDUCLAS2"] = "BKG" + header["HDUCLAS3"] = "FULL-ENCLOSURE" + header["HDUCLAS4"] = "BKG_2D" + header["FOVALIGN"] = "ALTAZ" + header["DATE"] = Time.now().utc.iso + _add_header_cards(header, **header_cards) + + return BinTableHDU(bkg, header=header, name=extname) + + @u.quantity_input( rad_max=u.deg, reco_energy_bins=u.TeV, diff --git a/pyirf/irf/__init__.py b/pyirf/irf/__init__.py index 3a7224684..82bae2372 100644 --- a/pyirf/irf/__init__.py +++ b/pyirf/irf/__init__.py @@ -2,16 +2,21 @@ 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 -from .background import background_2d +from .background import background_2d, background_3d_lonlat __all__ = [ "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", + "background_3d_lonlat", ] diff --git a/pyirf/irf/background.py b/pyirf/irf/background.py index 1e4561bd0..55b41f0ee 100644 --- a/pyirf/irf/background.py +++ b/pyirf/irf/background.py @@ -1,10 +1,10 @@ import astropy.units as u import numpy as np -from ..utils import cone_solid_angle +from ..utils import cone_solid_angle, rectangle_solid_angle #: Unit of the background rate IRF -BACKGROUND_UNIT = u.Unit('s-1 TeV-1 sr-1') +BACKGROUND_UNIT = u.Unit("s-1 TeV-1 sr-1") def background_2d(events, reco_energy_bins, fov_offset_bins, t_obs): @@ -43,7 +43,7 @@ def background_2d(events, reco_energy_bins, fov_offset_bins, t_obs): reco_energy_bins.to_value(u.TeV), fov_offset_bins.to_value(u.deg), ], - weights=events['weight'], + weights=events["weight"], ) # divide all energy bins by their width @@ -56,3 +56,69 @@ def background_2d(events, reco_energy_bins, fov_offset_bins, t_obs): bg_rate = per_energy / t_obs / bin_solid_angle return bg_rate.to(BACKGROUND_UNIT) + + +def background_3d_lonlat(events, reco_energy_bins, fov_lon_bins, fov_lat_bins, t_obs): + """ + Calculate background rates in square bins in the field of view. + + GADF documentation here: + https://gamma-astro-data-formats.readthedocs.io/en/latest/irfs/full_enclosure/bkg/index.html#bkg-3d + + Parameters + ---------- + events: astropy.table.QTable + DL2 events table of the selected background events. + Needed columns for this function: `reco_fov_lon`, `reco_fov_lat`, + `reco_energy`, `weight`. + reco_energy: astropy.units.Quantity[energy] + The bins in reconstructed energy to be used for the IRF + fov_lon_bins: astropy.units.Quantity[angle] + The bins in the field of view longitudinal direction to be used for the IRF. + Will become DETX bins. + fov_lat_bins: astropy.units.Quantity[angle] + The bins in the field of view latitudinal direction to be used for the IRF. + Will become DETY bins. + t_obs: astropy.units.Quantity[time] + Observation time. This must match with how the individual event + weights are calculated. + + Returns + ------- + bg_rate: astropy.units.Quantity + The background rate as particles per energy, time and solid angle + in the specified bins. + + Shape: (len(reco_energy_bins) - 1, len(fov_lon_bins) - 1, len(fov_lon_bins) - 1) + """ + fov_x_offset_bins = fov_lon_bins + fov_y_offset_bins = fov_lat_bins + + hist, _ = np.histogramdd( + [ + events["reco_energy"].to_value(u.TeV), + events["reco_fov_lon"].to_value(u.deg), + events["reco_fov_lat"].to_value(u.deg), + ], + bins=[ + reco_energy_bins.to_value(u.TeV), + fov_x_offset_bins.to_value(u.deg), + fov_y_offset_bins.to_value(u.deg), + ], + weights=events["weight"], + ) + + # divide all energy bins by their width + # hist has shape (n_energy, n_fov_offset) so we need to transpose and then back + bin_width_energy = np.diff(reco_energy_bins) + per_energy = (hist / bin_width_energy[:, np.newaxis, np.newaxis]) + + # divide by solid angle in each fov bin and the observation time + bin_solid_angle = rectangle_solid_angle( + fov_x_offset_bins[:-1], + fov_x_offset_bins[1:], + fov_y_offset_bins[:-1], + fov_y_offset_bins[1:], + ) + bg_rate = per_energy / t_obs / bin_solid_angle + return bg_rate.to(BACKGROUND_UNIT) 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/energy_dispersion.py b/pyirf/irf/energy_dispersion.py index de631e426..f34c3546d 100644 --- a/pyirf/irf/energy_dispersion.py +++ b/pyirf/irf/energy_dispersion.py @@ -6,7 +6,8 @@ __all__ = [ "energy_dispersion", - "energy_dispersion_to_migration" + "energy_migration_matrix", + "energy_dispersion_to_migration", ] @@ -29,7 +30,10 @@ def _normalize_hist(hist, migration_bins): def energy_dispersion( - selected_events, true_energy_bins, fov_offset_bins, migration_bins, + selected_events, + true_energy_bins, + fov_offset_bins, + migration_bins, ): """ Calculate energy dispersion for the given DL2 event list. @@ -80,6 +84,56 @@ def energy_dispersion( return energy_dispersion +@u.quantity_input(true_energy_bins=u.TeV, reco_energy_bins=u.TeV, fov_offset_bins=u.deg) +def energy_migration_matrix( + events, true_energy_bins, reco_energy_bins, fov_offset_bins +): + """Compute the energy migration matrix directly from the events. + + Parameters + ---------- + events : `~astropy.table.QTable` + Table of the DL2 events. + Required columns: ``reco_energy``, ``true_energy``, ``true_source_fov_offset``. + true_energy_bins : `~astropy.units.Quantity` + Bin edges in true energy. + reco_energy_bins : `~astropy.units.Quantity` + Bin edges in reconstructed energy. + + Returns + ------- + matrix : array-like + Migration matrix as probabilities along the reconstructed energy axis. + energy axis with shape + (n_true_energy_bins, n_reco_energy_bins, n_fov_offset_bins) + containing energies in TeV. + """ + + hist, _ = np.histogramdd( + np.column_stack( + [ + events["true_energy"].to_value(u.TeV), + events["reco_energy"].to_value(u.TeV), + events["true_source_fov_offset"].to_value(u.deg), + ] + ), + bins=[ + true_energy_bins.to_value(u.TeV), + reco_energy_bins.to_value(u.TeV), + fov_offset_bins.to_value(u.deg), + ], + ) + + with np.errstate(invalid="ignore"): + hist /= hist.sum(axis=1)[:, np.newaxis, :] + # the nans come from the fact that the sum along the reconstructed energy axis + # might sometimes be 0 when there are no events in that given true energy bin + # and fov offset bin + hist[np.isnan(hist)] = 0 + + return hist + + def energy_dispersion_to_migration( dispersion_matrix, disp_true_energy_edges, @@ -88,8 +142,8 @@ def energy_dispersion_to_migration( new_reco_energy_edges, ): """ - Construct a energy migration matrix from a energy - dispersion matrix. + Construct a energy migration matrix from an energy dispersion matrix. + Depending on the new energy ranges, the sum over the first axis can be smaller than 1. The new true energy bins need to be a subset of the old range, @@ -110,21 +164,25 @@ def energy_dispersion_to_migration( new_reco_energy_edges: astropy.units.Quantity[energy] Reco energy edges matching the second dimension of the output - Returns: - -------- + Returns + ------- migration_matrix: numpy.ndarray Three-dimensional energy migration matrix. The third dimension equals the fov offset dimension of the energy dispersion matrix. """ + migration_matrix = np.zeros( + ( + len(new_true_energy_edges) - 1, + len(new_reco_energy_edges) - 1, + dispersion_matrix.shape[2], + ) + ) - migration_matrix = np.zeros(( - len(new_true_energy_edges) - 1, - len(new_reco_energy_edges) - 1, - dispersion_matrix.shape[2], - )) + migra_width = np.diff(disp_migration_edges) + probability = dispersion_matrix * migra_width[np.newaxis, :, np.newaxis] true_energy_interpolation = resample_histogram1d( - dispersion_matrix, + probability, disp_true_energy_edges, new_true_energy_edges, axis=0, @@ -137,13 +195,14 @@ def energy_dispersion_to_migration( for idx, e_true in enumerate( (new_true_energy_edges[1:] + new_true_energy_edges[:-1]) / 2 ): - # get migration for the new true energy bin e_true_dispersion = true_energy_interpolation[idx] with warnings.catch_warnings(): # silence inf/inf division warning - warnings.filterwarnings('ignore', 'invalid value encountered in true_divide') + warnings.filterwarnings( + "ignore", "invalid value encountered in true_divide" + ) interpolation_edges = new_reco_energy_edges / e_true y = resample_histogram1d( diff --git a/pyirf/irf/tests/test_background.py b/pyirf/irf/tests/test_background.py index 23f760eb2..cfb4e2b18 100644 --- a/pyirf/irf/tests/test_background.py +++ b/pyirf/irf/tests/test_background.py @@ -36,4 +36,94 @@ def test_background(): # check that psf is normalized bin_solid_angle = np.diff(cone_solid_angle(fov_bins)) e_width = np.diff(energy_bins) - assert np.allclose(np.sum((bg.T * e_width).T * bin_solid_angle, axis=1), [1000, 100] / u.s) + assert np.allclose( + np.sum((bg.T * e_width).T * bin_solid_angle, axis=1), [1000, 100] / u.s + ) + + +def test_background_3d_lonlat(): + from pyirf.irf import background_3d_lonlat + from pyirf.utils import rectangle_solid_angle + from pyirf.irf.background import BACKGROUND_UNIT + + reco_energy_bins = [0.1, 1.1, 11.1, 111.1] * u.TeV + fov_lon_bins = [-1.0, 0, 1.0] * u.deg + fov_lat_bins = [-1.0, 0, 1.0] * u.deg + + N_low = 4000 + N_high = 40 + N_tot = N_low + N_high + + # Fill values + E_low, E_hig = 0.5, 5 + Lon_low, Lon_hig = (-0.5, 0.5) * u.deg + Lat_low, Lat_hig = (-0.5, 0.5) * u.deg + + t_obs = 100 * u.s + bin_width_energy = np.diff(reco_energy_bins) + bin_solid_angle = rectangle_solid_angle( + fov_lon_bins[:-1], fov_lon_bins[1:], fov_lat_bins[:-1], fov_lat_bins[1:] + ) + + # Toy events with two energies and four different sky positions + selected_events = QTable( + { + "reco_energy": np.concatenate( + [ + np.full(N_low // 4, E_low), + np.full(N_high // 4, E_hig), + np.full(N_low // 4, E_low), + np.full(N_high // 4, E_hig), + np.full(N_low // 4, E_low), + np.full(N_high // 4, E_hig), + np.full(N_low // 4, E_low), + np.full(N_high // 4, E_hig), + ] + ) + * u.TeV, + "reco_fov_lon": np.concatenate( + [ + np.full(N_low // 4, Lon_low), + np.full(N_high // 4, Lon_hig), + np.full(N_low // 4, Lon_low), + np.full(N_high // 4, Lon_hig), + np.full(N_low // 4, Lon_low), + np.full(N_high // 4, Lon_hig), + np.full(N_low // 4, Lon_low), + np.full(N_high // 4, Lon_hig), + ] + ) + * u.deg, + "reco_fov_lat": np.append( + np.full(N_tot // 2, Lat_low), + np.full(N_tot // 2, Lat_hig) + ) + * u.deg, + "weight": np.full(N_tot, 1.0), + } + ) + + bkg_rate = background_3d_lonlat( + selected_events, + reco_energy_bins=reco_energy_bins, + fov_lon_bins=fov_lon_bins, + fov_lat_bins=fov_lat_bins, + t_obs=t_obs, + ) + assert bkg_rate.shape == ( + len(reco_energy_bins) - 1, + len(fov_lon_bins) - 1, + len(fov_lat_bins) - 1, + ) + assert bkg_rate.unit == BACKGROUND_UNIT + + # Convert to counts, project to energy axis, and check counts round-trip correctly + assert np.allclose( + (bin_solid_angle * bkg_rate * bin_width_energy[:, np.newaxis, np.newaxis]).sum(axis=(1, 2)) * t_obs, + [N_low, N_high, 0], + ) + # Convert to counts, project to latitude axis, and check counts round-trip correctly + assert np.allclose( + (bin_solid_angle * bkg_rate * bin_width_energy[:, np.newaxis, np.newaxis]).sum(axis=(0, 1)) * t_obs, + 2 * [N_tot // 2], + ) 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/irf/tests/test_energy_dispersion.py b/pyirf/irf/tests/test_energy_dispersion.py index 036d106de..451511ffc 100644 --- a/pyirf/irf/tests/test_energy_dispersion.py +++ b/pyirf/irf/tests/test_energy_dispersion.py @@ -64,17 +64,29 @@ def ppf(cdf, bins, value): assert np.isclose( TRUE_SIGMA_1, - 0.5 * (ppf(cdf[0, :, 0], migration_bins, 0.84) - ppf(cdf[0, :, 0], migration_bins, 0.16)), + 0.5 + * ( + ppf(cdf[0, :, 0], migration_bins, 0.84) + - ppf(cdf[0, :, 0], migration_bins, 0.16) + ), rtol=0.1, ) assert np.isclose( TRUE_SIGMA_2, - 0.5 * (ppf(cdf[1, :, 0], migration_bins, 0.84) - ppf(cdf[1, :, 0], migration_bins, 0.16)), + 0.5 + * ( + ppf(cdf[1, :, 0], migration_bins, 0.84) + - ppf(cdf[1, :, 0], migration_bins, 0.16) + ), rtol=0.1, ) assert np.isclose( TRUE_SIGMA_3, - 0.5 * (ppf(cdf[2, :, 0], migration_bins, 0.84) - ppf(cdf[2, :, 0], migration_bins, 0.16)), + 0.5 + * ( + ppf(cdf[2, :, 0], migration_bins, 0.84) + - ppf(cdf[2, :, 0], migration_bins, 0.16) + ), rtol=0.1, ) @@ -85,16 +97,15 @@ def test_energy_dispersion_to_migration(): np.random.seed(0) N = 10000 - true_energy_bins = 10**np.arange(np.log10(0.2), np.log10(200), 1/10) * u.TeV + true_energy_bins = 10 ** np.arange(np.log10(0.2), np.log10(200), 1 / 10) * u.TeV fov_offset_bins = np.array([0, 1, 2]) * u.deg migration_bins = np.linspace(0, 2, 101) - true_energy = np.random.uniform( - true_energy_bins[0].value, - true_energy_bins[-1].value, - size=N - ) * u.TeV + true_energy = ( + np.random.uniform(true_energy_bins[0].value, true_energy_bins[-1].value, size=N) + * u.TeV + ) reco_energy = true_energy * np.random.uniform(0.5, 1.5, size=N) selected_events = QTable( @@ -116,15 +127,14 @@ def test_energy_dispersion_to_migration(): ) # migration matrix selecting a limited energy band with different binsizes - new_true_energy_bins = 10**np.arange(np.log10(2), np.log10(20), 1/5) * u.TeV - new_reco_energy_bins = 10**np.arange(np.log10(2), np.log10(20), 1/5) * u.TeV + new_true_energy_bins = 10 ** np.arange(np.log10(2), np.log10(20), 1 / 5) * u.TeV + new_reco_energy_bins = 10 ** np.arange(np.log10(2), np.log10(20), 1 / 5) * u.TeV migration_matrix = energy_dispersion_to_migration( dispersion_matrix, true_energy_bins, migration_bins, new_true_energy_bins, new_reco_energy_bins, - ) # test dimension @@ -137,10 +147,53 @@ def test_energy_dispersion_to_migration(): # test that migrations dont always sum to 1 (since some energies are # not included in the matrix) - assert migration_matrix.sum() < (len(new_true_energy_bins) - 1) * (len(fov_offset_bins) - 1) + assert migration_matrix.sum() < (len(new_true_energy_bins) - 1) * ( + len(fov_offset_bins) - 1 + ) assert np.all(np.isfinite(migration_matrix)) +def test_energy_migration_matrix_from_events(): + from pyirf.irf.energy_dispersion import energy_migration_matrix + + np.random.seed(0) + N = 10000 + true_energy_bins = 10 ** np.arange(np.log10(0.2), np.log10(200), 1 / 10) * u.TeV + reco_energy_bins = 10 ** np.arange(np.log10(2), np.log10(20), 1 / 5) * u.TeV + fov_offset_bins = np.array([0, 1, 2]) * u.deg + + true_energy = ( + np.random.uniform(true_energy_bins[0].value, true_energy_bins[-1].value, size=N) + * u.TeV + ) + reco_energy = true_energy * np.random.uniform(0.5, 1.5, size=N) + + events = QTable( + { + "reco_energy": reco_energy, + "true_energy": true_energy, + "true_source_fov_offset": np.concatenate( + [ + np.full(N // 2, 0.2), + np.full(N // 2, 1.5), + ] + ) + * u.deg, + } + ) + + matrix = energy_migration_matrix( + events, true_energy_bins, reco_energy_bins, fov_offset_bins + ) + + assert matrix.shape == ( + len(true_energy_bins) - 1, + len(reco_energy_bins) - 1, + len(fov_offset_bins) - 1, + ) + assert np.allclose(matrix.sum(axis=1).max(), 1, rtol=0.1) + + def test_overflow_bins_migration_matrix(): from pyirf.irf import energy_dispersion from pyirf.irf.energy_dispersion import energy_dispersion_to_migration @@ -150,21 +203,24 @@ def test_overflow_bins_migration_matrix(): N = 10000 # add under/overflow bins - bins = 10 ** np.arange( - np.log10(0.2), - np.log10(200), - 1 / 10, - ) * u.TeV + bins = ( + 10 + ** np.arange( + np.log10(0.2), + np.log10(200), + 1 / 10, + ) + * u.TeV + ) true_energy_bins = add_overflow_bins(bins, positive=False) fov_offset_bins = np.array([0, 1, 2]) * u.deg migration_bins = np.linspace(0, 2, 101) - true_energy = np.random.uniform( - true_energy_bins[1].value, - true_energy_bins[-2].value, - size=N - ) * u.TeV + true_energy = ( + np.random.uniform(true_energy_bins[1].value, true_energy_bins[-2].value, size=N) + * u.TeV + ) reco_energy = true_energy * np.random.uniform(0.5, 1.5, size=N) selected_events = QTable( diff --git a/pyirf/sensitivity.py b/pyirf/sensitivity.py index a8f7b114d..c8873c828 100644 --- a/pyirf/sensitivity.py +++ b/pyirf/sensitivity.py @@ -1,6 +1,7 @@ """ Functions to calculate sensitivity """ + import numpy as np from scipy.optimize import brentq import logging @@ -8,6 +9,7 @@ from astropy.table import QTable import astropy.units as u +from .compat import COPY_IF_NEEDED from .statistics import li_ma_significance from .utils import check_histograms, cone_solid_angle from .binning import create_histogram_table, bin_center @@ -36,7 +38,7 @@ def _relative_sensitivity( return np.nan if n_on < 0 or n_off < 0: - raise ValueError(f'n_on and n_off must be positive, got {n_on}, {n_off}') + raise ValueError(f"n_on and n_off must be positive, got {n_on}, {n_off}") n_background = n_off * alpha n_signal = n_on - n_background @@ -65,7 +67,7 @@ def equation(relative_flux): relative_flux = brentq(equation, lower_bound, upper_bound) except (RuntimeError, ValueError) as e: - log.warn( + log.warning( "Could not calculate relative significance for" f" n_signal={n_signal:.1f}, n_off={n_off:.1f}, returning nan {e}" ) @@ -87,8 +89,7 @@ def equation(relative_flux): _relative_sensitivity_vectorized = np.vectorize( - _relative_sensitivity, - excluded=['significance_function'] + _relative_sensitivity, excluded=["significance_function"] ) @@ -230,8 +231,10 @@ def calculate_sensitivity( # convert any quantities to arrays, # since quantitites don't work with vectorize - n_on = u.Quantity(n_on, copy=False).to_value(u.one) - n_off = u.Quantity(background_hist["n_weighted"], copy=False).to_value(u.one) + n_on = u.Quantity(n_on, copy=COPY_IF_NEEDED).to_value(u.one) + n_off = u.Quantity(background_hist["n_weighted"], copy=COPY_IF_NEEDED).to_value( + u.one + ) # calculate sensitivity in each bin rel_sens = relative_sensitivity( @@ -257,7 +260,9 @@ def calculate_sensitivity( s["n_background_weighted"] = background_hist["n_weighted"] # copy also "n_proton" / "n_electron_weighted" etc. if available - for k in filter(lambda c: c.startswith('n_') and c != 'n_weighted', background_hist.colnames): + for k in filter( + lambda c: c.startswith("n_") and c != "n_weighted", background_hist.colnames + ): s[k] = background_hist[k] s["significance"] = significance_function( @@ -271,10 +276,14 @@ def calculate_sensitivity( def estimate_background( - events, reco_energy_bins, theta_cuts, alpha, - fov_offset_min, fov_offset_max, + events, + reco_energy_bins, + theta_cuts, + alpha, + fov_offset_min, + fov_offset_max, ): - ''' + """ Estimate the number of background events for a point-like sensitivity. Due to limited statistics, it is often not possible to just apply the same @@ -308,33 +317,30 @@ def estimate_background( Minimum distance from the fov center for background events to be taken into account fov_offset_max: astropy.units.Quantity[angle] Maximum distance from the fov center for background events to be taken into account - ''' - in_ring = ( - (events['reco_source_fov_offset'] >= fov_offset_min) - & (events['reco_source_fov_offset'] < fov_offset_max) + """ + in_ring = (events["reco_source_fov_offset"] >= fov_offset_min) & ( + events["reco_source_fov_offset"] < fov_offset_max ) bg = create_histogram_table( events[in_ring], reco_energy_bins, - key='reco_energy', + key="reco_energy", ) # scale number of background events according to the on region size # background radius and alpha center = bin_center(reco_energy_bins) # interpolate the theta cut to the bins used here - theta_cuts_bg_bins = np.interp( - center, - theta_cuts['center'], - theta_cuts['cut'] - ) + theta_cuts_bg_bins = np.interp(center, theta_cuts["center"], theta_cuts["cut"]) - solid_angle_on = cone_solid_angle(theta_cuts_bg_bins) - solid_angle_ring = cone_solid_angle(fov_offset_max) - cone_solid_angle(fov_offset_min) + solid_angle_on = cone_solid_angle(theta_cuts_bg_bins) + solid_angle_ring = cone_solid_angle(fov_offset_max) - cone_solid_angle( + fov_offset_min + ) size_ratio = (solid_angle_on / solid_angle_ring).to_value(u.one) - for key in filter(lambda col: col.startswith('n'), bg.colnames): + for key in filter(lambda col: col.startswith("n"), bg.colnames): # *= not possible due to upcast from int to float bg[key] = bg[key] * size_ratio / alpha 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/spectral.py b/pyirf/spectral.py index 35bad8647..182657bff 100644 --- a/pyirf/spectral.py +++ b/pyirf/spectral.py @@ -1,6 +1,7 @@ """ Functions and classes for calculating spectral weights """ + import sys from importlib.resources import as_file, files @@ -9,6 +10,7 @@ from astropy.table import QTable from scipy.interpolate import interp1d +from .compat import COPY_IF_NEEDED from .utils import cone_solid_angle #: Unit of a point source flux @@ -120,7 +122,9 @@ def from_simulation(cls, simulated_event_info, obstime, e_ref=1 * u.TeV): viewcone_max = simulated_event_info.viewcone_max if (viewcone_max - viewcone_min).value > 0: - solid_angle = cone_solid_angle(viewcone_max) - cone_solid_angle(viewcone_min) + solid_angle = cone_solid_angle(viewcone_max) - cone_solid_angle( + viewcone_min + ) unit = DIFFUSE_FLUX_UNIT else: solid_angle = 1 @@ -308,7 +312,6 @@ def __init__( self.interp = interp1d(x, y, bounds_error=False, fill_value="extrapolate") def __call__(self, energy): - x = (energy / self.reference_energy).to_value(u.one) if self.log_energy: @@ -319,7 +322,7 @@ def __call__(self, energy): if self.log_flux: y = 10**y - return u.Quantity(y, self.flux_unit, copy=False) + return u.Quantity(y, self.flux_unit, copy=COPY_IF_NEEDED) @classmethod def from_table( diff --git a/pyirf/statistics.py b/pyirf/statistics.py index b4f5b10de..988bfd70b 100644 --- a/pyirf/statistics.py +++ b/pyirf/statistics.py @@ -1,5 +1,6 @@ import numpy as np +from .compat import COPY_IF_NEEDED from .utils import is_scalar __all__ = ["li_ma_significance"] @@ -34,8 +35,8 @@ def li_ma_significance(n_on, n_off, alpha=0.2): # Cast everything into float64 to avoid numeric instabilties # when multiplying very small and very big numbers to get t1 and t2 - n_on = np.array(n_on, copy=False, ndmin=1, dtype=np.float64) - n_off = np.array(n_off, copy=False, ndmin=1, dtype=np.float64) + n_on = np.array(n_on, copy=COPY_IF_NEEDED, ndmin=1, dtype=np.float64) + n_off = np.array(n_off, copy=COPY_IF_NEEDED, ndmin=1, dtype=np.float64) alpha = np.float64(alpha) with np.errstate(divide="ignore", invalid="ignore"): 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..e1d41c631 100644 --- a/pyirf/utils.py +++ b/pyirf/utils.py @@ -1,7 +1,9 @@ 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 .compat import COPY_IF_NEEDED from .exceptions import MissingColumns, WrongColumnUnit @@ -9,6 +11,7 @@ "is_scalar", "calculate_theta", "calculate_source_fov_offset", + "calculate_source_fov_position_angle", "check_histograms", "cone_solid_angle", ] @@ -27,7 +30,7 @@ def is_scalar(val): result: bool True is if input object is a scalar, False otherwise. """ - result = np.array(val, copy=False).shape == tuple() + result = np.array(val, copy=COPY_IF_NEEDED).shape == tuple() return result @@ -88,6 +91,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 +187,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. diff --git a/setup.py b/setup.py index 9c03055ab..e4ad60003 100644 --- a/setup.py +++ b/setup.py @@ -16,20 +16,17 @@ "notebook", "tables", "towncrier", - "astropy~=5.3", # gammapy doesn't yet support 6.0 but doesn't pin it gammapy, ], "tests": [ "pytest", "pytest-cov", - "astropy~=5.3", # gammapy doesn't yet support 6.0 but doesn't pin it gammapy, "ogadf-schema~=0.2.3", "uproot~=4.0", "awkward~=1.0", ], "gammapy": [ - "astropy~=5.3", # gammapy doesn't yet support 6.0 but doesn't pin it gammapy, ], } @@ -44,12 +41,13 @@ setup( use_scm_version={"write_to": os.path.join("pyirf", "_version.py")}, - packages=find_packages(exclude=['pyirf._dev_version']), + packages=find_packages(exclude=["pyirf._dev_version"]), install_requires=[ "astropy>=5.3,<7.0.0a0", "numpy>=1.21", - "scipy<1.12", + "scipy", "tqdm", + "packaging", ], include_package_data=True, extras_require=extras_require,