From 8042e5bb593365a8d828db867963aef13182ee1a Mon Sep 17 00:00:00 2001 From: Tomas Bylund Date: Tue, 16 Jan 2024 14:50:52 +0100 Subject: [PATCH 01/48] Add basic facilities for 3D background irfs --- pyirf/io/__init__.py | 2 ++ pyirf/io/gadf.py | 51 ++++++++++++++++++++++++++++++++++++++ pyirf/irf/__init__.py | 3 ++- pyirf/irf/background.py | 55 +++++++++++++++++++++++++++++++++++++++++ 4 files changed, 110 insertions(+), 1 deletion(-) 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 dd0d22f32..b34ff39c9 100644 --- a/pyirf/io/gadf.py +++ b/pyirf/io/gadf.py @@ -252,6 +252,57 @@ 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_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..37f2adfd0 100644 --- a/pyirf/irf/__init__.py +++ b/pyirf/irf/__init__.py @@ -5,7 +5,7 @@ ) from .energy_dispersion import energy_dispersion from .psf import psf_table -from .background import background_2d +from .background import background_2d, background_3d __all__ = [ "effective_area", @@ -14,4 +14,5 @@ "energy_dispersion", "psf_table", "background_2d", + "background_3d", ] diff --git a/pyirf/irf/background.py b/pyirf/irf/background.py index 1e4561bd0..2d050ca84 100644 --- a/pyirf/irf/background.py +++ b/pyirf/irf/background.py @@ -56,3 +56,58 @@ 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(events, reco_energy_bins, fov_offset_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_offset_bins: astropy.units.Quantity[angle] + The bins in the field of view offset to be used for the IRF + 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_offset_bins) - 1, len(fov_offset_bins) - 1) + """ + 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_offset_bins.to_value(u.deg), + fov_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.T / bin_width_energy).T + + # divide by solid angle in each fov bin and the observation time + bin_solid_angle = np.diff(fov_offset_bins) + bg_rate = per_energy / t_obs / bin_solid_angle**2 + + return bg_rate.to(BACKGROUND_UNIT) From 7c6d15848131e29ea7405596e5ef29481d00b151 Mon Sep 17 00:00:00 2001 From: Tomas Bylund Date: Tue, 6 Feb 2024 15:26:27 +0100 Subject: [PATCH 02/48] Should now support rectangular arrays --- pyirf/irf/background.py | 21 ++++++++++++++++----- 1 file changed, 16 insertions(+), 5 deletions(-) diff --git a/pyirf/irf/background.py b/pyirf/irf/background.py index 2d050ca84..1c4242d12 100644 --- a/pyirf/irf/background.py +++ b/pyirf/irf/background.py @@ -4,7 +4,7 @@ from ..utils import cone_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 @@ -74,7 +74,7 @@ def background_3d(events, reco_energy_bins, fov_offset_bins, t_obs): reco_energy: astropy.units.Quantity[energy] The bins in reconstructed energy to be used for the IRF fov_offset_bins: astropy.units.Quantity[angle] - The bins in the field of view offset to be used for the IRF + The bins in the field of view offset to be used for the IRF, either a (N,) or (1,N) or a (2,N) array t_obs: astropy.units.Quantity[time] Observation time. This must match with how the individual event weights are calculated. @@ -87,6 +87,17 @@ def background_3d(events, reco_energy_bins, fov_offset_bins, t_obs): Shape: (len(reco_energy_bins) - 1, len(fov_offset_bins) - 1, len(fov_offset_bins) - 1) """ + if (fov_offset_bins.shape[0] == 1) or (len(fov_offset_bins.shape) == 1): + fov_x_offset_bins = fov_offset_bins + fov_y_offset_bins = fov_offset_bins + elif fov_offset_bins.shape[0] == 2: + fov_x_offset_bins = fov_offset_bins[0, :] + fov_y_offset_bins = fov_offset_bins[1, :] + else: + raise ValueError( + f"fov_offset_bins must be eiher (N,) or (1,N) or (2,N), found {fov_offset_bins.shape}" + ) + hist, _ = np.histogramdd( [ events["reco_energy"].to_value(u.TeV), @@ -95,8 +106,8 @@ def background_3d(events, reco_energy_bins, fov_offset_bins, t_obs): ], bins=[ reco_energy_bins.to_value(u.TeV), - fov_offset_bins.to_value(u.deg), - fov_offset_bins.to_value(u.deg), + fov_x_offset_bins.to_value(u.deg), + fov_y_offset_bins.to_value(u.deg), ], weights=events["weight"], ) From 2990605798e5f10c03dbac0ca9a37694c04f8639 Mon Sep 17 00:00:00 2001 From: Maximilian Linhoff Date: Tue, 2 Apr 2024 16:38:44 +0200 Subject: [PATCH 03/48] Remove compatibility pins --- setup.py | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/setup.py b/setup.py index 9c03055ab..561acc80d 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, ], } @@ -48,7 +45,7 @@ install_requires=[ "astropy>=5.3,<7.0.0a0", "numpy>=1.21", - "scipy<1.12", + "scipy", "tqdm", ], include_package_data=True, From 17a75148ebd6d9b6c8073c9aac2879afc448d5b7 Mon Sep 17 00:00:00 2001 From: Rune Michael Dominik Date: Mon, 8 Apr 2024 16:07:18 +0200 Subject: [PATCH 04/48] Add 2023 ICRC paper to README --- README.rst | 17 +++++++++++++++++ 1 file changed, 17 insertions(+) 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}, + } From ec666904451281360ba3a744d13a93d570dad426 Mon Sep 17 00:00:00 2001 From: Tomas Bylund Date: Tue, 16 Jan 2024 14:50:52 +0100 Subject: [PATCH 05/48] Add basic facilities for 3D background irfs --- pyirf/io/__init__.py | 2 ++ pyirf/io/gadf.py | 51 ++++++++++++++++++++++++++++++++++++++ pyirf/irf/__init__.py | 3 ++- pyirf/irf/background.py | 55 +++++++++++++++++++++++++++++++++++++++++ 4 files changed, 110 insertions(+), 1 deletion(-) 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..824d00839 100644 --- a/pyirf/io/gadf.py +++ b/pyirf/io/gadf.py @@ -260,6 +260,57 @@ 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_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..37f2adfd0 100644 --- a/pyirf/irf/__init__.py +++ b/pyirf/irf/__init__.py @@ -5,7 +5,7 @@ ) from .energy_dispersion import energy_dispersion from .psf import psf_table -from .background import background_2d +from .background import background_2d, background_3d __all__ = [ "effective_area", @@ -14,4 +14,5 @@ "energy_dispersion", "psf_table", "background_2d", + "background_3d", ] diff --git a/pyirf/irf/background.py b/pyirf/irf/background.py index 1e4561bd0..2d050ca84 100644 --- a/pyirf/irf/background.py +++ b/pyirf/irf/background.py @@ -56,3 +56,58 @@ 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(events, reco_energy_bins, fov_offset_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_offset_bins: astropy.units.Quantity[angle] + The bins in the field of view offset to be used for the IRF + 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_offset_bins) - 1, len(fov_offset_bins) - 1) + """ + 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_offset_bins.to_value(u.deg), + fov_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.T / bin_width_energy).T + + # divide by solid angle in each fov bin and the observation time + bin_solid_angle = np.diff(fov_offset_bins) + bg_rate = per_energy / t_obs / bin_solid_angle**2 + + return bg_rate.to(BACKGROUND_UNIT) From 6d4119738e34fd4b74f09603838e52eb72b7ee95 Mon Sep 17 00:00:00 2001 From: Tomas Bylund Date: Tue, 6 Feb 2024 15:26:27 +0100 Subject: [PATCH 06/48] Should now support rectangular arrays --- pyirf/irf/background.py | 21 ++++++++++++++++----- 1 file changed, 16 insertions(+), 5 deletions(-) diff --git a/pyirf/irf/background.py b/pyirf/irf/background.py index 2d050ca84..1c4242d12 100644 --- a/pyirf/irf/background.py +++ b/pyirf/irf/background.py @@ -4,7 +4,7 @@ from ..utils import cone_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 @@ -74,7 +74,7 @@ def background_3d(events, reco_energy_bins, fov_offset_bins, t_obs): reco_energy: astropy.units.Quantity[energy] The bins in reconstructed energy to be used for the IRF fov_offset_bins: astropy.units.Quantity[angle] - The bins in the field of view offset to be used for the IRF + The bins in the field of view offset to be used for the IRF, either a (N,) or (1,N) or a (2,N) array t_obs: astropy.units.Quantity[time] Observation time. This must match with how the individual event weights are calculated. @@ -87,6 +87,17 @@ def background_3d(events, reco_energy_bins, fov_offset_bins, t_obs): Shape: (len(reco_energy_bins) - 1, len(fov_offset_bins) - 1, len(fov_offset_bins) - 1) """ + if (fov_offset_bins.shape[0] == 1) or (len(fov_offset_bins.shape) == 1): + fov_x_offset_bins = fov_offset_bins + fov_y_offset_bins = fov_offset_bins + elif fov_offset_bins.shape[0] == 2: + fov_x_offset_bins = fov_offset_bins[0, :] + fov_y_offset_bins = fov_offset_bins[1, :] + else: + raise ValueError( + f"fov_offset_bins must be eiher (N,) or (1,N) or (2,N), found {fov_offset_bins.shape}" + ) + hist, _ = np.histogramdd( [ events["reco_energy"].to_value(u.TeV), @@ -95,8 +106,8 @@ def background_3d(events, reco_energy_bins, fov_offset_bins, t_obs): ], bins=[ reco_energy_bins.to_value(u.TeV), - fov_offset_bins.to_value(u.deg), - fov_offset_bins.to_value(u.deg), + fov_x_offset_bins.to_value(u.deg), + fov_y_offset_bins.to_value(u.deg), ], weights=events["weight"], ) From 7fc4c741e8114419ee696a33d811b733b752d4c4 Mon Sep 17 00:00:00 2001 From: Tomas Bylund Date: Fri, 8 Nov 2024 15:20:22 +0100 Subject: [PATCH 07/48] Fixes requested by Max, added test --- pyirf/io/gadf.py | 19 +++++-- pyirf/irf/background.py | 10 ++-- pyirf/irf/tests/test_background.py | 91 +++++++++++++++++++++++++++++- 3 files changed, 109 insertions(+), 11 deletions(-) diff --git a/pyirf/io/gadf.py b/pyirf/io/gadf.py index 824d00839..963b6de78 100644 --- a/pyirf/io/gadf.py +++ b/pyirf/io/gadf.py @@ -266,7 +266,9 @@ def create_background_2d_hdu( def create_background_3d_hdu( background_3d, reco_energy_bins, - fov_offset_bins, + fov_lon_bins, + fov_lat_bins, + alignment = "ALTAZ", extname="BACKGROUND", **header_cards, ): @@ -282,8 +284,13 @@ def create_background_3d_hdu( (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. + 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 + Wheter 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 @@ -293,8 +300,8 @@ def create_background_3d_hdu( 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)) + 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) @@ -304,7 +311,7 @@ def create_background_3d_hdu( header["HDUCLAS2"] = "BKG" header["HDUCLAS3"] = "FULL-ENCLOSURE" header["HDUCLAS4"] = "BKG_2D" - header["FOVALIGN"] = "ALTAZ" + header["FOVALIGN"] = alignment header["DATE"] = Time.now().utc.iso _add_header_cards(header, **header_cards) diff --git a/pyirf/irf/background.py b/pyirf/irf/background.py index 1c4242d12..a5f9b3b0d 100644 --- a/pyirf/irf/background.py +++ b/pyirf/irf/background.py @@ -1,7 +1,7 @@ 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") @@ -118,7 +118,9 @@ def background_3d(events, reco_energy_bins, fov_offset_bins, t_obs): per_energy = (hist.T / bin_width_energy).T # divide by solid angle in each fov bin and the observation time - bin_solid_angle = np.diff(fov_offset_bins) - bg_rate = per_energy / t_obs / bin_solid_angle**2 - + 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/tests/test_background.py b/pyirf/irf/tests/test_background.py index 23f760eb2..5d3e009ae 100644 --- a/pyirf/irf/tests/test_background.py +++ b/pyirf/irf/tests/test_background.py @@ -36,4 +36,93 @@ 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(): + from pyirf.irf import background_3d + 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( + selected_events, + reco_energy_bins=reco_energy_bins, + fov_offset_bins=np.vstack((fov_lat_bins, fov_lon_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.T * bin_width_energy).T).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.T * bin_width_energy).T).sum(axis=(0, 1)) * t_obs, + 2 * [N_tot // 2], + ) From a3325f3cba8713e62f6d852d85f81a7cc2de0f38 Mon Sep 17 00:00:00 2001 From: Rune Michael Dominik Date: Fri, 16 Feb 2024 14:09:12 +0100 Subject: [PATCH 08/48] Add tests for more correct handling of fill-values in RadMaxEstimator --- ...st_component_estimator_specific_classes.py | 191 +++++++++++++++++- 1 file changed, 189 insertions(+), 2 deletions(-) diff --git a/pyirf/interpolation/tests/test_component_estimator_specific_classes.py b/pyirf/interpolation/tests/test_component_estimator_specific_classes.py index c1f77913b..dbca00b51 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,185 @@ 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]]) + + estim = 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, + ) + + with pytest.raises( + ValueError, + match="Fill-value handling only supported in up to two grid dimensions.", + ): + estim(np.array([0.25, 0.25, 0.25])) From 0b3e1779fe317da76dcc4a626b8592bce77052ae Mon Sep 17 00:00:00 2001 From: Rune Michael Dominik Date: Fri, 16 Feb 2024 14:10:01 +0100 Subject: [PATCH 09/48] Add more correct handling of fill-values in RadMaxEstimator --- pyirf/interpolation/component_estimators.py | 120 +++++++++++++++++++- 1 file changed, 118 insertions(+), 2 deletions(-) diff --git a/pyirf/interpolation/component_estimators.py b/pyirf/interpolation/component_estimators.py index 740180e60..f45b6c1e2 100644 --- a/pyirf/interpolation/component_estimators.py +++ b/pyirf/interpolation/component_estimators.py @@ -3,7 +3,6 @@ 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 @@ -13,7 +12,9 @@ PDFNormalization, ) from .griddata_interpolator import GridDataInterpolator +from .nearest_neighbor_searcher import BaseNearestNeighborSearcher from .quantile_interpolator import QuantileInterpolator +from .utils import find_nearest_simplex __all__ = [ "BaseComponentEstimator", @@ -477,6 +478,7 @@ def __init__( self, grid_points, rad_max, + fill_value=None, interpolator_cls=GridDataInterpolator, interpolator_kwargs=None, extrapolator_cls=None, @@ -495,6 +497,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 +532,20 @@ 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 + + # 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 +563,101 @@ 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 + full_estimation[full_estimation < 0] = 0 + + # Early exit if fill_values should not be handled + if not self.fill_val: + return full_estimation + + # Raise error if fill-values should be handled in >=3 dims + if self.grid_dim >= 3: + raise ValueError( + "Fill-value handling only supported in up to two grid dimensions." + ) + + # Early exit if a nearest neighbor estimation would be overwritten + 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): From 7d8f79f8484c319201a898463ae594caff2f536d Mon Sep 17 00:00:00 2001 From: Rune Michael Dominik Date: Thu, 21 Mar 2024 13:43:35 +0100 Subject: [PATCH 10/48] Add newsfragment --- docs/changes/282.feature.rst | 1 + 1 file changed, 1 insertion(+) create mode 100644 docs/changes/282.feature.rst diff --git a/docs/changes/282.feature.rst b/docs/changes/282.feature.rst new file mode 100644 index 000000000..3c61c6a44 --- /dev/null +++ b/docs/changes/282.feature.rst @@ -0,0 +1 @@ +Add basic combinatoric fill-value handling for RAD_MAX estimation. From 211f7128a3dba685fc3c34c85550dd97d0ed66fa Mon Sep 17 00:00:00 2001 From: Rune Michael Dominik Date: Mon, 8 Apr 2024 12:00:14 +0200 Subject: [PATCH 11/48] adress comments --- pyirf/interpolation/component_estimators.py | 20 +++++++++++-------- ...st_component_estimator_specific_classes.py | 20 +++++++++---------- 2 files changed, 21 insertions(+), 19 deletions(-) diff --git a/pyirf/interpolation/component_estimators.py b/pyirf/interpolation/component_estimators.py index f45b6c1e2..cd14beade 100644 --- a/pyirf/interpolation/component_estimators.py +++ b/pyirf/interpolation/component_estimators.py @@ -541,6 +541,12 @@ def __init__( 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): @@ -566,19 +572,15 @@ def __call__(self, target_point): # First, construct estimation without handling fill-values full_estimation = super().__call__(target_point) # Safeguard against extreme extrapolation cases - full_estimation[full_estimation < 0] = 0 + 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 - # Raise error if fill-values should be handled in >=3 dims - if self.grid_dim >= 3: - raise ValueError( - "Fill-value handling only supported in up to two grid dimensions." - ) - # 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()) @@ -631,7 +633,9 @@ def __call__(self, target_point): # 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") + 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 diff --git a/pyirf/interpolation/tests/test_component_estimator_specific_classes.py b/pyirf/interpolation/tests/test_component_estimator_specific_classes.py index dbca00b51..af786c4e6 100644 --- a/pyirf/interpolation/tests/test_component_estimator_specific_classes.py +++ b/pyirf/interpolation/tests/test_component_estimator_specific_classes.py @@ -358,18 +358,16 @@ def test_RadMaxEstimator_fill_val_handling_3D(): rad_max = np.array([[0.95], [0.95], [0.95], [0.95]]) - estim = 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, - ) - with pytest.raises( ValueError, match="Fill-value handling only supported in up to two grid dimensions.", ): - estim(np.array([0.25, 0.25, 0.25])) + 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, + ) From a46b73a0cbb0ce146db44a8796f2f2edfe1d60ee Mon Sep 17 00:00:00 2001 From: Maximilian Linhoff Date: Tue, 14 May 2024 12:06:20 +0200 Subject: [PATCH 12/48] Render changelog for 0.11 --- CHANGES.rst | 37 ++++++++++++++++++++++++++++++++ docs/changes/253.feature.rst | 2 -- docs/changes/264.feature.rst | 1 - docs/changes/266.maintenance.rst | 1 - docs/changes/268.bugfix.rst | 4 ---- docs/changes/271.maintenance.rst | 1 - docs/changes/275.maintenance.rst | 1 - docs/changes/279.maintenance.rst | 1 - docs/changes/282.feature.rst | 1 - 9 files changed, 37 insertions(+), 12 deletions(-) delete mode 100644 docs/changes/253.feature.rst delete mode 100644 docs/changes/264.feature.rst delete mode 100644 docs/changes/266.maintenance.rst delete mode 100644 docs/changes/268.bugfix.rst delete mode 100644 docs/changes/271.maintenance.rst delete mode 100644 docs/changes/275.maintenance.rst delete mode 100644 docs/changes/279.maintenance.rst delete mode 100644 docs/changes/282.feature.rst 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/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/268.bugfix.rst b/docs/changes/268.bugfix.rst deleted file mode 100644 index 7458041a1..000000000 --- a/docs/changes/268.bugfix.rst +++ /dev/null @@ -1,4 +0,0 @@ -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. 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/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/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/282.feature.rst b/docs/changes/282.feature.rst deleted file mode 100644 index 3c61c6a44..000000000 --- a/docs/changes/282.feature.rst +++ /dev/null @@ -1 +0,0 @@ -Add basic combinatoric fill-value handling for RAD_MAX estimation. From 16829149b526736b2942518ee78ab333f9a2f976 Mon Sep 17 00:00:00 2001 From: Lukas Beiske Date: Mon, 30 Sep 2024 15:43:59 +0200 Subject: [PATCH 13/48] Make things compatible with numpy 2.0; replace deprecated logging.warn() --- pyirf/benchmarks/energy_bias_resolution.py | 11 +++-- pyirf/binning.py | 36 ++++++++------ pyirf/compat.py | 10 ++++ pyirf/cuts.py | 26 +++++++---- pyirf/interpolation/component_estimators.py | 15 ++++-- pyirf/sensitivity.py | 52 ++++++++++++--------- pyirf/spectral.py | 9 ++-- pyirf/statistics.py | 5 +- pyirf/utils.py | 3 +- setup.py | 3 +- 10 files changed, 107 insertions(+), 63 deletions(-) create mode 100644 pyirf/compat.py 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/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/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 cd14beade..06e123f37 100644 --- a/pyirf/interpolation/component_estimators.py +++ b/pyirf/interpolation/component_estimators.py @@ -1,16 +1,19 @@ """Classes to estimate (interpolate/extrapolate) actual IRF HDUs""" + import warnings import astropy.units as u import numpy as np 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 @@ -428,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__( @@ -468,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): @@ -852,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/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/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/utils.py b/pyirf/utils.py index d5ed730e5..28e1da5eb 100644 --- a/pyirf/utils.py +++ b/pyirf/utils.py @@ -2,6 +2,7 @@ import astropy.units as u from astropy.coordinates import angular_separation +from .compat import COPY_IF_NEEDED from .exceptions import MissingColumns, WrongColumnUnit @@ -27,7 +28,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 diff --git a/setup.py b/setup.py index 561acc80d..e4ad60003 100644 --- a/setup.py +++ b/setup.py @@ -41,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", "tqdm", + "packaging", ], include_package_data=True, extras_require=extras_require, From dd68f0616d4ee6f24dff3d9612249f59c413a598 Mon Sep 17 00:00:00 2001 From: Lukas Beiske Date: Mon, 30 Sep 2024 15:53:05 +0200 Subject: [PATCH 14/48] Add changelog entry --- docs/changes/289.maintenance.rst | 3 +++ 1 file changed, 3 insertions(+) create mode 100644 docs/changes/289.maintenance.rst 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``. From fa5af1f0b0acbba9d8d70b2f572ff9553700e010 Mon Sep 17 00:00:00 2001 From: Michele Peresano Date: Tue, 28 May 2024 13:54:35 +0200 Subject: [PATCH 15/48] Update .mailmap --- .mailmap | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) 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 From 0803c77ed756cfd7cb3aa393cc8caf4ac9526bf3 Mon Sep 17 00:00:00 2001 From: Lukas Beiske Date: Mon, 30 Sep 2024 17:51:59 +0200 Subject: [PATCH 16/48] Add option of multiple quantiles to angular_resolution --- docs/changes/290.feature.rst | 2 + pyirf/benchmarks/angular_resolution.py | 30 ++++++--- .../tests/test_angular_resolution.py | 67 ++++++++++++------- 3 files changed, 68 insertions(+), 31 deletions(-) create mode 100644 docs/changes/290.feature.rst diff --git a/docs/changes/290.feature.rst b/docs/changes/290.feature.rst new file mode 100644 index 000000000..42dacd1b5 --- /dev/null +++ b/docs/changes/290.feature.rst @@ -0,0 +1,2 @@ +Make it possible to pass multiple quantiles to ``pyirf.benchmarks.angular_resolution``, calculating all of them. +If only one quantile is passed, the output is the same as before. diff --git a/pyirf/benchmarks/angular_resolution.py b/pyirf/benchmarks/angular_resolution.py index 2f3445cfe..fe07e3c47 100644 --- a/pyirf/benchmarks/angular_resolution.py +++ b/pyirf/benchmarks/angular_resolution.py @@ -10,9 +10,10 @@ def angular_resolution( - events, energy_bins, + events, + energy_bins, energy_type="true", - quantile=ONE_SIGMA_QUANTILE, + quantile=[ONE_SIGMA_QUANTILE], ): """ Calculate the angular resolution. @@ -20,6 +21,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 +32,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 +55,17 @@ 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) + # keep backwards compatibility + if not isinstance(quantile, list) and not isinstance(quantile, np.ndarray): + quantile = [quantile] + + if len(quantile) < 2: + keys = ["angular_resolution"] + else: + keys = [f"containment_quantile_{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 +75,8 @@ 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) + for key, value in zip(keys, quantile): + result[key][bin_idx] = np.nanquantile(group["theta"], value) + result["n_events"][bin_idx] = len(group) + return result diff --git a/pyirf/benchmarks/tests/test_angular_resolution.py b/pyirf/benchmarks/tests/test_angular_resolution.py index 269c6c46b..fc4941ae8 100644 --- a/pyirf/benchmarks/tests/test_angular_resolution.py +++ b/pyirf/benchmarks/tests/test_angular_resolution.py @@ -8,18 +8,18 @@ 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, + } ) + table = angular_resolution(events, [1, 10, 100] * u.TeV) + assert np.all(np.isnan(table["angular_resolution"])) + @pytest.mark.parametrize("unit", (u.deg, u.rad)) def test_angular_resolution(unit): from pyirf.benchmarks import angular_resolution @@ -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"].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"].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["containment_quantile_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["containment_quantile_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["containment_quantile_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) From 254c4c19bd27cc976ce1cb3f55ee33a2cea1eef6 Mon Sep 17 00:00:00 2001 From: Lukas Beiske Date: Tue, 1 Oct 2024 13:51:40 +0200 Subject: [PATCH 17/48] No mutable default argument; check for sequence, not list and np.ndarray --- pyirf/benchmarks/angular_resolution.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/pyirf/benchmarks/angular_resolution.py b/pyirf/benchmarks/angular_resolution.py index fe07e3c47..577b9e588 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 @@ -13,7 +14,7 @@ def angular_resolution( events, energy_bins, energy_type="true", - quantile=[ONE_SIGMA_QUANTILE], + quantile=ONE_SIGMA_QUANTILE, ): """ Calculate the angular resolution. @@ -55,8 +56,7 @@ def angular_resolution( result[f"{energy_key}_center"] = 0.5 * (energy_bins[:-1] + energy_bins[1:]) result["n_events"] = 0 - # keep backwards compatibility - if not isinstance(quantile, list) and not isinstance(quantile, np.ndarray): + if not isinstance(quantile, Sequence): quantile = [quantile] if len(quantile) < 2: From 592c4c6462a3021dd6ea1b0a71129fd1f0fa043c Mon Sep 17 00:00:00 2001 From: Lukas Beiske Date: Tue, 8 Oct 2024 14:25:53 +0200 Subject: [PATCH 18/48] Do not change column naming based on number of quantiles --- examples/comparison_with_EventDisplay.ipynb | 2 +- pyirf/benchmarks/angular_resolution.py | 5 +---- pyirf/benchmarks/tests/test_angular_resolution.py | 6 +++--- 3 files changed, 5 insertions(+), 8 deletions(-) diff --git a/examples/comparison_with_EventDisplay.ipynb b/examples/comparison_with_EventDisplay.ipynb index 71c2d59c4..aa6321e28 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['containment_quantile_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 577b9e588..0fde76a9d 100644 --- a/pyirf/benchmarks/angular_resolution.py +++ b/pyirf/benchmarks/angular_resolution.py @@ -59,10 +59,7 @@ def angular_resolution( if not isinstance(quantile, Sequence): quantile = [quantile] - if len(quantile) < 2: - keys = ["angular_resolution"] - else: - keys = [f"containment_quantile_{value * 100:.0f}" for value in quantile] + keys = [f"containment_quantile_{value * 100:.0f}" for value in quantile] for key in keys: result[key] = u.Quantity(np.nan, table["theta"].unit) diff --git a/pyirf/benchmarks/tests/test_angular_resolution.py b/pyirf/benchmarks/tests/test_angular_resolution.py index fc4941ae8..07a7ce661 100644 --- a/pyirf/benchmarks/tests/test_angular_resolution.py +++ b/pyirf/benchmarks/tests/test_angular_resolution.py @@ -17,7 +17,7 @@ def test_empty_angular_resolution(): table = angular_resolution(events, [1, 10, 100] * u.TeV) - assert np.all(np.isnan(table["angular_resolution"])) + assert np.all(np.isnan(table["containment_quantile_68"])) @pytest.mark.parametrize("unit", (u.deg, u.rad)) @@ -57,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["containment_quantile_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) @@ -68,7 +68,7 @@ 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["containment_quantile_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) From 38fff116eb6ec74ab42bbccdeb8dd47fb706d6e7 Mon Sep 17 00:00:00 2001 From: Maximilian Linhoff Date: Mon, 21 Oct 2024 14:59:04 +0200 Subject: [PATCH 19/48] Fix docs warning --- docs/conf.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) 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 = { From b80d8fb52c51607a036fe007b35f563a49f7cbc5 Mon Sep 17 00:00:00 2001 From: Lukas Beiske Date: Mon, 21 Oct 2024 15:59:07 +0200 Subject: [PATCH 20/48] Address comments --- docs/changes/290.feature.rst | 1 - examples/comparison_with_EventDisplay.ipynb | 2 +- pyirf/benchmarks/angular_resolution.py | 2 +- pyirf/benchmarks/tests/test_angular_resolution.py | 12 ++++++------ 4 files changed, 8 insertions(+), 9 deletions(-) diff --git a/docs/changes/290.feature.rst b/docs/changes/290.feature.rst index 42dacd1b5..2be73b37e 100644 --- a/docs/changes/290.feature.rst +++ b/docs/changes/290.feature.rst @@ -1,2 +1 @@ Make it possible to pass multiple quantiles to ``pyirf.benchmarks.angular_resolution``, calculating all of them. -If only one quantile is passed, the output is the same as before. diff --git a/examples/comparison_with_EventDisplay.ipynb b/examples/comparison_with_EventDisplay.ipynb index aa6321e28..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['containment_quantile_68'].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 0fde76a9d..a8d4f146e 100644 --- a/pyirf/benchmarks/angular_resolution.py +++ b/pyirf/benchmarks/angular_resolution.py @@ -59,7 +59,7 @@ def angular_resolution( if not isinstance(quantile, Sequence): quantile = [quantile] - keys = [f"containment_quantile_{value * 100:.0f}" for value in 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) diff --git a/pyirf/benchmarks/tests/test_angular_resolution.py b/pyirf/benchmarks/tests/test_angular_resolution.py index 07a7ce661..e66ed50dd 100644 --- a/pyirf/benchmarks/tests/test_angular_resolution.py +++ b/pyirf/benchmarks/tests/test_angular_resolution.py @@ -17,7 +17,7 @@ def test_empty_angular_resolution(): table = angular_resolution(events, [1, 10, 100] * u.TeV) - assert np.all(np.isnan(table["containment_quantile_68"])) + assert np.all(np.isnan(table["angular_resolution_68"])) @pytest.mark.parametrize("unit", (u.deg, u.rad)) @@ -57,7 +57,7 @@ def test_angular_resolution(unit): bins = [1, 10, 100] * u.TeV table = angular_resolution(events, bins) - ang_res = table["containment_quantile_68"].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) @@ -68,14 +68,14 @@ 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["containment_quantile_95"].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["containment_quantile_25"].to(u.deg) + 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 @@ -84,10 +84,10 @@ def test_angular_resolution(unit): cov_25[1], norm(0, TRUE_RES_2).interval(0.25)[1] * u.deg, rtol=0.05 ) - cov_50 = table["containment_quantile_50"].to(u.deg) + 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["containment_quantile_90"].to(u.deg) + 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) From fe9e4787bd15f24978843a9dfe0c8640034a3ea2 Mon Sep 17 00:00:00 2001 From: Lukas Beiske Date: Mon, 21 Oct 2024 18:02:40 +0200 Subject: [PATCH 21/48] Call np.nanquantile with multiple quantiles --- pyirf/benchmarks/angular_resolution.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/pyirf/benchmarks/angular_resolution.py b/pyirf/benchmarks/angular_resolution.py index a8d4f146e..d87f60ba9 100644 --- a/pyirf/benchmarks/angular_resolution.py +++ b/pyirf/benchmarks/angular_resolution.py @@ -72,8 +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): - for key, value in zip(keys, quantile): - result[key][bin_idx] = np.nanquantile(group["theta"], value) - result["n_events"][bin_idx] = len(group) + 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 From 01fb541c4952e3fbbece807a58dc91708c270b7f Mon Sep 17 00:00:00 2001 From: Lukas Beiske Date: Wed, 6 Nov 2024 11:05:00 +0100 Subject: [PATCH 22/48] Rename and update changelog --- docs/changes/290.api.rst | 3 +++ docs/changes/290.feature.rst | 1 - 2 files changed, 3 insertions(+), 1 deletion(-) create mode 100644 docs/changes/290.api.rst delete mode 100644 docs/changes/290.feature.rst 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/changes/290.feature.rst b/docs/changes/290.feature.rst deleted file mode 100644 index 2be73b37e..000000000 --- a/docs/changes/290.feature.rst +++ /dev/null @@ -1 +0,0 @@ -Make it possible to pass multiple quantiles to ``pyirf.benchmarks.angular_resolution``, calculating all of them. From 84fda49c6f24f87349bd63f1508743918be14fe7 Mon Sep 17 00:00:00 2001 From: Luca Di Bella Date: Mon, 12 Feb 2024 17:35:58 +0100 Subject: [PATCH 23/48] Add 2 methods to calculate effective area in energy and 2 spacial dimensions The different methods use binning in different coordinate systems - effective_area_3D_polar uses a polar coordinate system with radial FOV offset and azimuthal FOV position angle bins - effective_area_3D_nominal uses a nominal coordinate system with FOV longitude and FOV latitude bins --- pyirf/irf/effective_area.py | 66 +++++++++++++++++++++++++++++++++++++ 1 file changed, 66 insertions(+) diff --git a/pyirf/irf/effective_area.py b/pyirf/irf/effective_area.py index daa13e767..6ad542f28 100644 --- a/pyirf/irf/effective_area.py +++ b/pyirf/irf/effective_area.py @@ -85,3 +85,69 @@ def effective_area_per_energy_and_fov( ) return effective_area(hist_selected, hist_simulated, area) + +def effective_area_3D_polar( + selected_events, simulation_info, energy_bins, fov_offset_bins, fov_position_angle_bins +): + """ + Calculate effective area in bins of true energy, field of view offset, and field of view position angle. + + Parameters + ---------- + selected_events: astropy.table.QTable + DL2 events table, required columns for this function: + - `true_energy` + - `true_source_fov_offset` + - `true_source_fov_position_angle` + simulation_info: pyirf.simulations.SimulatedEventsInfo + The overall statistics of the simulated events + true_energy_bins: astropy.units.Quantity[energy] + The true energy bin edges in which to calculate effective area. + fov_offset_bins: astropy.units.Quantity[angle] + The field of view radial bin edges in which to calculate effective area. + fov_position_angle_bins: astropy.units.Quantity[radian] + The field of view azimuthal bin edges in which to calculate effective area. + """ + area = np.pi * simulation_info.max_impact ** 2 + + hist_simulated = simulation_info.calculate_n_showers_3D_radial(simulation_info, energy_bins, fov_offset_bins, fov_position_angle_bins) + + hist_selected,_ = np.histogramdd( + np.array([selected_events['true_energy'].to_value(u.TeV), selected_events['true_source_fov_offset'].to_value(u.deg), selected_events['true_source_fov_position_angle'].to_value(u.rad)]).T, + bins=(energy_bins.to_value(u.TeV), fov_offset_bins.to_value(u.deg), fov_position_angle_bins.to_value(u.rad)), + ) + + return effective_area(hist_selected, hist_simulated, area) + +def effective_area_3D_nominal( + selected_events, simulation_info, energy_bins, fov_longitude_bins, fov_latitude_bins +): + """ + Calculate effective area in bins of true energy, field of view longitude, and field of view latitude. + + Parameters + ---------- + selected_events: astropy.table.QTable + DL2 events table, required columns for this function: + - `true_energy` + - `true_source_fov_lon` + - `true_source_fov_lat` + simulation_info: pyirf.simulations.SimulatedEventsInfo + The overall statistics of the simulated events + true_energy_bins: astropy.units.Quantity[energy] + The true energy bin edges in which to calculate effective area. + fov_longitude_bins: astropy.units.Quantity[angle] + The field of view longitude bin edges in which to calculate effective area. + fov_latitude_bins: astropy.units.Quantity[angle] + The field of view latitude bin edges in which to calculate effective area. + """ + area = np.pi * simulation_info.max_impact ** 2 + + hist_simulated = simulation_info.calculate_n_showers_3D_nominal(simulation_info, energy_bins, fov_longitude_bins, fov_latitude_bins) + + hist_selected, _ = np.histogramdd( + np.array([selected_events['true_energy'].to_value(u.TeV), selected_events['true_source_fov_lon'].value, selected_events['true_source_fov_lat'].value]).T, + bins=(energy_bins.to_value(u.TeV), fov_longitude_bins.to_value(u.deg), fov_latitude_bins.to_value(u.deg)), + ) + + return effective_area(hist_selected, hist_simulated, area) From 68c3fa635347c24e081ae8eda8be47b4a27d7149 Mon Sep 17 00:00:00 2001 From: Luca Di Bella Date: Mon, 12 Feb 2024 17:46:38 +0100 Subject: [PATCH 24/48] Add methods to calculate n_showers in energy and 2 spacial dimensions These calculate the number of showers in each bin for different coordinate systems - calculate_n_showers_3D_polar uses a polar coordinate system with radial FOV offset bins and azimuthal FOV position angle bins - calculate_n_showers_3D_nominal uses a nominal coordinate system with FOV longitude and latitude bins --- pyirf/simulations.py | 76 ++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 76 insertions(+) diff --git a/pyirf/simulations.py b/pyirf/simulations.py index 2bedf4c1b..b3b1b2773 100644 --- a/pyirf/simulations.py +++ b/pyirf/simulations.py @@ -1,5 +1,6 @@ import astropy.units as u import numpy as np +from .utils import cone_solid_angle __all__ = [ 'SimulatedEventsInfo', @@ -168,6 +169,81 @@ def calculate_n_showers_per_energy_and_fov(self, energy_bins, fov_bins): fov_integral = self.calculate_n_showers_per_fov(fov_bins) return e_integral[:, np.newaxis] * fov_integral / self.n_showers + @u.quantity_input(energy_bins=u.TeV, fov_offset_bins=u.deg, fov_position_angle_bins=u.rad) + def calculate_n_showers_3D_polar(self, energy_bins, fov_offset_bins, fov_position_angle_bins): + """ + Calculate number of showers that were simulated in the given + energy and 2D fov bins in polar coordinates. + + This assumes the events were generated uniformly distributed per solid angle, + and from a powerlaw in energy like CORSIKA simulates events. + + Parameters + ---------- + energy_bins: astropy.units.Quantity[energy] + The energy bin edges for which to calculate the number of simulated showers + fov_offset_bins: astropy.units.Quantity[angle] + The FOV radial bin edges for which to calculate the number of simulated showers + fov_position_angle_bins: astropy.units.Quantity[radian] + The FOV azimuthal bin edges for which to calculate the number of simulated showers + + + Returns + ------- + n_showers: numpy.ndarray(ndim=3) + The expected number of events inside each of the + ``energy_bins``, ``fov_offset_bins`` and ``fov_position_angle_bins``. + Dimension (n_energy_bins, n_fov_offset_bins, n_fov_position_angle_bins) + This is a floating point number. + The actual numbers will follow a poissionian distribution around this + expected value. + """ + e_fov_offset_integral = self.calculate_n_showers_per_energy_and_fov(energy_bins, fov_offset_bins) + position_angle_integral = np.ones(len(fov_position_angle_bins) - 1) * self.calculate_n_showers_per_fov(np.linspace(self.viewcone_min, self.viewcone_max, 2)) / (len(fov_position_angle_bins) - 1) + + return e_fov_offset_integral[:,:,np.newaxis] * position_angle_integral / self.n_showers + + @u.quantity_input(energy_bins=u.TeV, fov_longitude_bins=u.deg, fov_latitude_bins=u.rad) + def calculate_n_showers_3D_nominal(self, energy_bins, fov_longitude_bins, fov_latitude_bins): + """ + Calculate number of showers that were simulated in the given + energy and 2D fov bins in nominal coordinates. + + This assumes the events were generated uniformly distributed per solid angle, + and from a powerlaw in energy like CORSIKA simulates events. + + Parameters + ---------- + energy_bins: astropy.units.Quantity[energy] + The energy bin edges for which to calculate the number of simulated showers + fov_longitude_bins: astropy.units.Quantity[angle] + The FOV longitude bin edges for which to calculate the number of simulated showers + fov_latitude_bins: astropy.units.Quantity[angle] + The FOV latitude bin edges for which to calculate the number of simulated showers + + + Returns + ------- + n_showers: numpy.ndarray(ndim=3) + The expected number of events inside each of the + ``energy_bins``, ``fov_longitude_bins`` and ``fov_latitude_bins``. + Dimension (n_energy_bins, n_fov_longitude_bins, n_fov_latitude_bins) + This is a floating point number. + The actual numbers will follow a poissionian distribution around this + expected value. + """ + fov_bin_midpoints_lon, fov_bin_midpoints_lat = np.meshgrid((fov_longitude_bins[:-1]+fov_longitude_bins[1:])/2,(fov_latitude_bins[:-1]+fov_latitude_bins[1:])/2) + mask_outside_fov = np.logical_or( np.sqrt( fov_bin_midpoints_lon**2 + fov_bin_midpoints_lat**2 ) > self.viewcone_max, np.sqrt( fov_bin_midpoints_lon**2 + fov_bin_midpoints_lat**2 ) < self.viewcone_min) + A_bins = np.outer(np.diff( np.sin( fov_latitude_bins.to_value(u.rad) ) ), np.diff(fov_longitude_bins.to_value(u.rad))) * u.sr + shower_density = self.n_showers / (cone_solid_angle(self.viewcone_max) - cone_solid_angle(self.viewcone_min)) + + fov_integral = shower_density * A_bins + e_integral = self.calculate_n_showers_per_energy(energy_bins) + + fov_integral[mask_outside_fov] = 0 + + return (e_integral[:,np.newaxis,np.newaxis] * fov_integral) / self.n_showers + def __repr__(self): return ( f"{self.__class__.__name__}(" From 242ecab25686b54ec6b905e567c6c2ab122837b8 Mon Sep 17 00:00:00 2001 From: Luca Di Bella Date: Tue, 13 Feb 2024 15:31:00 +0100 Subject: [PATCH 25/48] Add util for calculating FOV position angle --- pyirf/utils.py | 29 +++++++++++++++++++++++++++++ 1 file changed, 29 insertions(+) diff --git a/pyirf/utils.py b/pyirf/utils.py index 28e1da5eb..22dab61c6 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 astropy.coordinates import position_angle from .compat import COPY_IF_NEEDED from .exceptions import MissingColumns, WrongColumnUnit @@ -89,6 +90,34 @@ def calculate_source_fov_offset(events, prefix="true"): return theta.to(u.deg) +def calculate_source_fov_position_angle(events, prefix="true"): + """Calculate position_angle of true positions relative to pointing positions. + + Parameters + ---------- + events : astropy.QTable + Astropy Table object containing the reconstructed events information. + + prefix: str + Column prefix for az / alt, can be used to calculate reco or true + source fov position_angle. + + Returns + ------- + phi: astropy.units.Quantity + Position angle of the true positions relative to the pointing positions + in the sky. + """ + phi = position_angle( + events["pointing_az"], + events["pointing_alt"], + events[f"{prefix}_az"], + events[f"{prefix}_alt"], + ) + + return phi.to(u.deg) + + def check_histograms(hist1, hist2, key="reco_energy"): """ Check if two histogram tables have the same binning From 50d8002a2712c8d9815ccf10a210f0f07706029d Mon Sep 17 00:00:00 2001 From: Luca Di Bella Date: Tue, 13 Feb 2024 15:53:34 +0100 Subject: [PATCH 26/48] Fix effective_area_3d and poition angle util --- pyirf/irf/__init__.py | 4 ++++ pyirf/irf/effective_area.py | 6 ++++-- pyirf/simulations.py | 4 ++-- pyirf/utils.py | 1 + 4 files changed, 11 insertions(+), 4 deletions(-) diff --git a/pyirf/irf/__init__.py b/pyirf/irf/__init__.py index 37f2adfd0..5d5cb08cc 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_nominal, ) from .energy_dispersion import energy_dispersion from .psf import psf_table @@ -11,6 +13,8 @@ "effective_area", "effective_area_per_energy", "effective_area_per_energy_and_fov", + "effective_area_3d_polar", + "effective_area_3d_nominal", "energy_dispersion", "psf_table", "background_2d", diff --git a/pyirf/irf/effective_area.py b/pyirf/irf/effective_area.py index 6ad542f28..452e0f6cf 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_nominal", ] @@ -86,7 +88,7 @@ def effective_area_per_energy_and_fov( return effective_area(hist_selected, hist_simulated, area) -def effective_area_3D_polar( +def effective_area_3d_polar( selected_events, simulation_info, energy_bins, fov_offset_bins, fov_position_angle_bins ): """ @@ -119,7 +121,7 @@ def effective_area_3D_polar( return effective_area(hist_selected, hist_simulated, area) -def effective_area_3D_nominal( +def effective_area_3d_nominal( selected_events, simulation_info, energy_bins, fov_longitude_bins, fov_latitude_bins ): """ diff --git a/pyirf/simulations.py b/pyirf/simulations.py index b3b1b2773..e327801e5 100644 --- a/pyirf/simulations.py +++ b/pyirf/simulations.py @@ -170,7 +170,7 @@ def calculate_n_showers_per_energy_and_fov(self, energy_bins, 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): + 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. @@ -204,7 +204,7 @@ def calculate_n_showers_3D_polar(self, energy_bins, fov_offset_bins, fov_positio return e_fov_offset_integral[:,:,np.newaxis] * position_angle_integral / self.n_showers @u.quantity_input(energy_bins=u.TeV, fov_longitude_bins=u.deg, fov_latitude_bins=u.rad) - def calculate_n_showers_3D_nominal(self, energy_bins, fov_longitude_bins, fov_latitude_bins): + def calculate_n_showers_3d_nominal(self, energy_bins, fov_longitude_bins, fov_latitude_bins): """ Calculate number of showers that were simulated in the given energy and 2D fov bins in nominal coordinates. diff --git a/pyirf/utils.py b/pyirf/utils.py index 22dab61c6..d5ee765a5 100644 --- a/pyirf/utils.py +++ b/pyirf/utils.py @@ -11,6 +11,7 @@ "is_scalar", "calculate_theta", "calculate_source_fov_offset", + "calculate_source_fov_position_angle" "check_histograms", "cone_solid_angle", ] From 6547e22ce611853060b74ec501b504eaf46a104e Mon Sep 17 00:00:00 2001 From: Luca Di Bella Date: Tue, 13 Feb 2024 16:00:36 +0100 Subject: [PATCH 27/48] Fix missing comma --- pyirf/utils.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyirf/utils.py b/pyirf/utils.py index d5ee765a5..4be243521 100644 --- a/pyirf/utils.py +++ b/pyirf/utils.py @@ -11,7 +11,7 @@ "is_scalar", "calculate_theta", "calculate_source_fov_offset", - "calculate_source_fov_position_angle" + "calculate_source_fov_position_angle", "check_histograms", "cone_solid_angle", ] From b43c0aa1bc0898c6b8a9b2b255f4be725624e1f0 Mon Sep 17 00:00:00 2001 From: Luca Di Bella Date: Thu, 15 Feb 2024 13:54:41 +0100 Subject: [PATCH 28/48] Fix typos causing error when calling calculate_n_showers_3d_* --- pyirf/irf/effective_area.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pyirf/irf/effective_area.py b/pyirf/irf/effective_area.py index 452e0f6cf..20beeb507 100644 --- a/pyirf/irf/effective_area.py +++ b/pyirf/irf/effective_area.py @@ -112,7 +112,7 @@ def effective_area_3d_polar( """ area = np.pi * simulation_info.max_impact ** 2 - hist_simulated = simulation_info.calculate_n_showers_3D_radial(simulation_info, energy_bins, fov_offset_bins, fov_position_angle_bins) + hist_simulated = simulation_info.calculate_n_showers_3d_polar(simulation_info, energy_bins, fov_offset_bins, fov_position_angle_bins) hist_selected,_ = np.histogramdd( np.array([selected_events['true_energy'].to_value(u.TeV), selected_events['true_source_fov_offset'].to_value(u.deg), selected_events['true_source_fov_position_angle'].to_value(u.rad)]).T, @@ -145,7 +145,7 @@ def effective_area_3d_nominal( """ area = np.pi * simulation_info.max_impact ** 2 - hist_simulated = simulation_info.calculate_n_showers_3D_nominal(simulation_info, energy_bins, fov_longitude_bins, fov_latitude_bins) + hist_simulated = simulation_info.calculate_n_showers_3d_nominal(simulation_info, energy_bins, fov_longitude_bins, fov_latitude_bins) hist_selected, _ = np.histogramdd( np.array([selected_events['true_energy'].to_value(u.TeV), selected_events['true_source_fov_lon'].value, selected_events['true_source_fov_lat'].value]).T, From a7c5a9a0cd293a11db8efff9e4774cb5dd435b56 Mon Sep 17 00:00:00 2001 From: Luca Di Bella Date: Thu, 15 Feb 2024 14:00:43 +0100 Subject: [PATCH 29/48] Remove redundant argument breaking calculate_n_showers_3d call --- pyirf/irf/effective_area.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pyirf/irf/effective_area.py b/pyirf/irf/effective_area.py index 20beeb507..4f45d0081 100644 --- a/pyirf/irf/effective_area.py +++ b/pyirf/irf/effective_area.py @@ -112,7 +112,7 @@ def effective_area_3d_polar( """ area = np.pi * simulation_info.max_impact ** 2 - hist_simulated = simulation_info.calculate_n_showers_3d_polar(simulation_info, energy_bins, fov_offset_bins, fov_position_angle_bins) + hist_simulated = simulation_info.calculate_n_showers_3d_polar(energy_bins, fov_offset_bins, fov_position_angle_bins) hist_selected,_ = np.histogramdd( np.array([selected_events['true_energy'].to_value(u.TeV), selected_events['true_source_fov_offset'].to_value(u.deg), selected_events['true_source_fov_position_angle'].to_value(u.rad)]).T, @@ -145,7 +145,7 @@ def effective_area_3d_nominal( """ area = np.pi * simulation_info.max_impact ** 2 - hist_simulated = simulation_info.calculate_n_showers_3d_nominal(simulation_info, energy_bins, fov_longitude_bins, fov_latitude_bins) + hist_simulated = simulation_info.calculate_n_showers_3d_nominal(energy_bins, fov_longitude_bins, fov_latitude_bins) hist_selected, _ = np.histogramdd( np.array([selected_events['true_energy'].to_value(u.TeV), selected_events['true_source_fov_lon'].value, selected_events['true_source_fov_lat'].value]).T, From e033135d62c26b93abf7d32a3079cf0ad253c5fe Mon Sep 17 00:00:00 2001 From: Luca Di Bella Date: Thu, 15 Feb 2024 17:13:06 +0100 Subject: [PATCH 30/48] Fix formatting and issues from pull request Effective Area 3D (no. 281) code review --- pyirf/irf/__init__.py | 4 +-- pyirf/irf/effective_area.py | 66 +++++++++++++++++++++++++---------- pyirf/simulations.py | 68 +++++++++++++++++++++++++++---------- 3 files changed, 101 insertions(+), 37 deletions(-) diff --git a/pyirf/irf/__init__.py b/pyirf/irf/__init__.py index 5d5cb08cc..b002b4289 100644 --- a/pyirf/irf/__init__.py +++ b/pyirf/irf/__init__.py @@ -3,7 +3,7 @@ effective_area_per_energy_and_fov, effective_area_per_energy, effective_area_3d_polar, - effective_area_3d_nominal, + effective_area_3d_lonlat, ) from .energy_dispersion import energy_dispersion from .psf import psf_table @@ -14,7 +14,7 @@ "effective_area_per_energy", "effective_area_per_energy_and_fov", "effective_area_3d_polar", - "effective_area_3d_nominal", + "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 4f45d0081..102a3add6 100644 --- a/pyirf/irf/effective_area.py +++ b/pyirf/irf/effective_area.py @@ -8,7 +8,7 @@ "effective_area_per_energy", "effective_area_per_energy_and_fov", "effective_area_3d_polar", - "effective_area_3d_nominal", + "effective_area_3d_lonlat", ] @@ -88,8 +88,13 @@ def effective_area_per_energy_and_fov( return effective_area(hist_selected, hist_simulated, area) + def effective_area_3d_polar( - selected_events, simulation_info, energy_bins, fov_offset_bins, fov_position_angle_bins + selected_events, + simulation_info, + energy_bins, + fov_offset_bins, + fov_position_angle_bins, ): """ Calculate effective area in bins of true energy, field of view offset, and field of view position angle. @@ -110,18 +115,31 @@ def effective_area_3d_polar( fov_position_angle_bins: astropy.units.Quantity[radian] The field of view azimuthal bin edges in which to calculate effective area. """ - area = np.pi * simulation_info.max_impact ** 2 - - hist_simulated = simulation_info.calculate_n_showers_3d_polar(energy_bins, fov_offset_bins, fov_position_angle_bins) - - hist_selected,_ = np.histogramdd( - np.array([selected_events['true_energy'].to_value(u.TeV), selected_events['true_source_fov_offset'].to_value(u.deg), selected_events['true_source_fov_position_angle'].to_value(u.rad)]).T, - bins=(energy_bins.to_value(u.TeV), fov_offset_bins.to_value(u.deg), fov_position_angle_bins.to_value(u.rad)), - ) - + area = np.pi * simulation_info.max_impact**2 + + hist_simulated = simulation_info.calculate_n_showers_3d_polar( + energy_bins, fov_offset_bins, fov_position_angle_bins + ) + + hist_selected, _ = np.histogramdd( + np.column_stack( + [ + selected_events["true_energy"].to_value(u.TeV), + selected_events["true_source_fov_offset"].to_value(u.deg), + selected_events["true_source_fov_position_angle"].to_value(u.rad), + ] + ), + bins=( + energy_bins.to_value(u.TeV), + fov_offset_bins.to_value(u.deg), + fov_position_angle_bins.to_value(u.rad), + ), + ) + return effective_area(hist_selected, hist_simulated, area) -def effective_area_3d_nominal( + +def effective_area_3d_lonlat( selected_events, simulation_info, energy_bins, fov_longitude_bins, fov_latitude_bins ): """ @@ -143,13 +161,25 @@ def effective_area_3d_nominal( fov_latitude_bins: astropy.units.Quantity[angle] The field of view latitude bin edges in which to calculate effective area. """ - area = np.pi * simulation_info.max_impact ** 2 - - hist_simulated = simulation_info.calculate_n_showers_3d_nominal(energy_bins, fov_longitude_bins, fov_latitude_bins) - + area = np.pi * simulation_info.max_impact**2 + + hist_simulated = simulation_info.calculate_n_showers_3d_lonlat( + energy_bins, fov_longitude_bins, fov_latitude_bins + ) + hist_selected, _ = np.histogramdd( - np.array([selected_events['true_energy'].to_value(u.TeV), selected_events['true_source_fov_lon'].value, selected_events['true_source_fov_lat'].value]).T, - bins=(energy_bins.to_value(u.TeV), fov_longitude_bins.to_value(u.deg), fov_latitude_bins.to_value(u.deg)), + np.column_stack( + [ + selected_events["true_energy"].to_value(u.TeV), + selected_events["true_source_fov_lon"].to_value(u.deg), + selected_events["true_source_fov_lat"].to_value(u.deg), + ] + ), + bins=( + energy_bins.to_value(u.TeV), + fov_longitude_bins.to_value(u.deg), + fov_latitude_bins.to_value(u.deg), + ), ) return effective_area(hist_selected, hist_simulated, area) diff --git a/pyirf/simulations.py b/pyirf/simulations.py index e327801e5..89913087d 100644 --- a/pyirf/simulations.py +++ b/pyirf/simulations.py @@ -169,8 +169,12 @@ 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): + @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. @@ -187,7 +191,6 @@ def calculate_n_showers_3d_polar(self, energy_bins, fov_offset_bins, fov_positio 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) @@ -198,13 +201,29 @@ def calculate_n_showers_3d_polar(self, energy_bins, fov_offset_bins, fov_positio The actual numbers will follow a poissionian distribution around this expected value. """ - e_fov_offset_integral = self.calculate_n_showers_per_energy_and_fov(energy_bins, fov_offset_bins) - position_angle_integral = np.ones(len(fov_position_angle_bins) - 1) * self.calculate_n_showers_per_fov(np.linspace(self.viewcone_min, self.viewcone_max, 2)) / (len(fov_position_angle_bins) - 1) + e_fov_offset_integral = self.calculate_n_showers_per_energy_and_fov( + energy_bins, fov_offset_bins + ) + position_angle_integral = np.full( + len(fov_position_angle_bins) - 1, + self.calculate_n_showers_per_fov( + np.linspace(self.viewcone_min, self.viewcone_max, 2) + ) + / (len(fov_position_angle_bins) - 1), + ) - return e_fov_offset_integral[:,:,np.newaxis] * position_angle_integral / self.n_showers - - @u.quantity_input(energy_bins=u.TeV, fov_longitude_bins=u.deg, fov_latitude_bins=u.rad) - def calculate_n_showers_3d_nominal(self, energy_bins, fov_longitude_bins, fov_latitude_bins): + return ( + e_fov_offset_integral[:, :, np.newaxis] + * position_angle_integral + / self.n_showers + ) + + @u.quantity_input( + energy_bins=u.TeV, fov_longitude_bins=u.deg, fov_latitude_bins=u.rad + ) + def calculate_n_showers_3d_lonlat( + self, energy_bins, fov_longitude_bins, fov_latitude_bins + ): """ Calculate number of showers that were simulated in the given energy and 2D fov bins in nominal coordinates. @@ -221,7 +240,6 @@ def calculate_n_showers_3d_nominal(self, energy_bins, fov_longitude_bins, fov_la 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) @@ -232,17 +250,33 @@ def calculate_n_showers_3d_nominal(self, energy_bins, fov_longitude_bins, fov_la The actual numbers will follow a poissionian distribution around this expected value. """ - fov_bin_midpoints_lon, fov_bin_midpoints_lat = np.meshgrid((fov_longitude_bins[:-1]+fov_longitude_bins[1:])/2,(fov_latitude_bins[:-1]+fov_latitude_bins[1:])/2) - mask_outside_fov = np.logical_or( np.sqrt( fov_bin_midpoints_lon**2 + fov_bin_midpoints_lat**2 ) > self.viewcone_max, np.sqrt( fov_bin_midpoints_lon**2 + fov_bin_midpoints_lat**2 ) < self.viewcone_min) - A_bins = np.outer(np.diff( np.sin( fov_latitude_bins.to_value(u.rad) ) ), np.diff(fov_longitude_bins.to_value(u.rad))) * u.sr - shower_density = self.n_showers / (cone_solid_angle(self.viewcone_max) - cone_solid_angle(self.viewcone_min)) - + fov_bin_midpoints_lon, fov_bin_midpoints_lat = np.meshgrid( + (fov_longitude_bins[:-1] + fov_longitude_bins[1:]) / 2, + (fov_latitude_bins[:-1] + fov_latitude_bins[1:]) / 2, + ) + mask_outside_fov = np.logical_or( + np.sqrt(fov_bin_midpoints_lon**2 + fov_bin_midpoints_lat**2) + > self.viewcone_max, + np.sqrt(fov_bin_midpoints_lon**2 + fov_bin_midpoints_lat**2) + < self.viewcone_min, + ) + A_bins = ( + np.outer( + np.diff(np.sin(fov_latitude_bins.to_value(u.rad))), + np.diff(fov_longitude_bins.to_value(u.rad)), + ) + * u.sr + ) + shower_density = self.n_showers / ( + cone_solid_angle(self.viewcone_max) - cone_solid_angle(self.viewcone_min) + ) + fov_integral = shower_density * A_bins e_integral = self.calculate_n_showers_per_energy(energy_bins) fov_integral[mask_outside_fov] = 0 - - return (e_integral[:,np.newaxis,np.newaxis] * fov_integral) / self.n_showers + + return (e_integral[:, np.newaxis, np.newaxis] * fov_integral) / self.n_showers def __repr__(self): return ( From cb58e414cd938b1ebceaf5421a85d4d617baed34 Mon Sep 17 00:00:00 2001 From: Luca Di Bella Date: Fri, 1 Mar 2024 15:58:57 +0100 Subject: [PATCH 31/48] Add functions to transform from sky coordinates to FOV coordinates --- pyirf/coordinates.py | 62 +++++++++++++++++++++++++++++++++ pyirf/tests/test_coordinates.py | 60 +++++++++++++++++++++++++++++++ 2 files changed, 122 insertions(+) create mode 100644 pyirf/coordinates.py create mode 100644 pyirf/tests/test_coordinates.py diff --git a/pyirf/coordinates.py b/pyirf/coordinates.py new file mode 100644 index 000000000..96be28054 --- /dev/null +++ b/pyirf/coordinates.py @@ -0,0 +1,62 @@ +import astropy.units as u +from astropy.coordinates import SkyCoord, SkyOffsetFrame, angular_separation, position_angle + +__all__ = [ + 'fov_coords_lon_lat', + 'fov_coords_theta_phi' +] + + +def fov_coords_lon_lat(lon, lat, pointing_lon, pointing_lat): + """Transform sky coordinates to field-of-view longitude-latitude coordinates. + + 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 fov_coords_theta_phi(az, alt, pointing_az, pointing_alt): + """Transform sky coordinates to field-of-view theta-phi coordinates. + + Parameters + ---------- + az, alt : `~astropy.units.Quantity` + Sky coordinate to be transformed. + pointing_az, pointing_alt : `~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. + """ + + theta = angular_separation(pointing_az, pointing_alt, az, alt) + + phi = position_angle(pointing_az, pointing_alt, az, alt) + + return theta.to(u.deg), (-phi).wrap_at(360 * u.deg).to(u.deg) \ 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..5a04a5d19 --- /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_fov_coords_lon_lat(): + from pyirf.coordinates import fov_coords_lon_lat + # test some simple cases + lon, lat = 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 = 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 = 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 = 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_fov_coords_theta_phi(): + from pyirf.coordinates import fov_coords_theta_phi + + theta, phi = fov_coords_theta_phi( + alt=1 * u.deg, az=0 * u.deg, pointing_alt=0 * u.deg, pointing_az=0 * u.deg + ) + assert u.isclose(theta, 1 * u.deg) + assert u.isclose(phi, 0 * u.deg) + + theta, phi = fov_coords_theta_phi( + alt=-1 * u.deg, az=0 * u.deg, pointing_alt=0 * u.deg, pointing_az=0 * u.deg + ) + assert u.isclose(theta, 1 * u.deg) + assert u.isclose(phi, 180 * u.deg) + + theta, phi = fov_coords_theta_phi( + alt=0 * u.deg, az=-1 * u.deg, pointing_alt=0 * u.deg, pointing_az=0 * u.deg + ) + assert u.isclose(theta, 1 * u.deg) + assert u.isclose(phi, 90 * u.deg) + + theta, phi = fov_coords_theta_phi( + alt=0 * u.deg, az=1 * u.deg, pointing_alt=0 * u.deg, pointing_az=0 * u.deg + ) + assert u.isclose(theta, 1 * u.deg) + assert u.isclose(phi, 270 * u.deg) \ No newline at end of file From 2ec0a19752d523b456be804ae4441eefe92aeaa6 Mon Sep 17 00:00:00 2001 From: Luca Di Bella Date: Fri, 1 Mar 2024 16:47:24 +0100 Subject: [PATCH 32/48] Update calculate_n_showers_3d functions and add tests --- pyirf/simulations.py | 125 ++++++++++++++++++++++++-------- pyirf/tests/test_simulations.py | 88 ++++++++++++++++++++++ 2 files changed, 181 insertions(+), 32 deletions(-) diff --git a/pyirf/simulations.py b/pyirf/simulations.py index 89913087d..fc99bdc78 100644 --- a/pyirf/simulations.py +++ b/pyirf/simulations.py @@ -1,6 +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', @@ -204,25 +207,22 @@ def calculate_n_showers_3d_polar( e_fov_offset_integral = self.calculate_n_showers_per_energy_and_fov( energy_bins, fov_offset_bins ) - position_angle_integral = np.full( - len(fov_position_angle_bins) - 1, - self.calculate_n_showers_per_fov( - np.linspace(self.viewcone_min, self.viewcone_max, 2) - ) - / (len(fov_position_angle_bins) - 1), + 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 ( - e_fov_offset_integral[:, :, np.newaxis] - * position_angle_integral - / self.n_showers - ) + 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 + self, energy_bins, fov_longitude_bins, fov_latitude_bins, subpixels=20 ): """ Calculate number of showers that were simulated in the given @@ -250,31 +250,26 @@ def calculate_n_showers_3d_lonlat( The actual numbers will follow a poissionian distribution around this expected value. """ - fov_bin_midpoints_lon, fov_bin_midpoints_lat = np.meshgrid( - (fov_longitude_bins[:-1] + fov_longitude_bins[1:]) / 2, - (fov_latitude_bins[:-1] + fov_latitude_bins[1:]) / 2, - ) - mask_outside_fov = np.logical_or( - np.sqrt(fov_bin_midpoints_lon**2 + fov_bin_midpoints_lat**2) - > self.viewcone_max, - np.sqrt(fov_bin_midpoints_lon**2 + fov_bin_midpoints_lat**2) - < self.viewcone_min, - ) - A_bins = ( - np.outer( - np.diff(np.sin(fov_latitude_bins.to_value(u.rad))), - np.diff(fov_longitude_bins.to_value(u.rad)), - ) - * u.sr + + fov_mask = _fov_lonlat_grid_overlap_mask( + self, fov_longitude_bins, fov_latitude_bins, self.viewcone_max, subpixels=subpixels ) - shower_density = self.n_showers / ( - cone_solid_angle(self.viewcone_max) - cone_solid_angle(self.viewcone_min) + + 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 * A_bins + fov_integral = shower_density * bin_area e_integral = self.calculate_n_showers_per_energy(energy_bins) - fov_integral[mask_outside_fov] = 0 + fov_integral = fov_mask * fov_integral return (e_integral[:, np.newaxis, np.newaxis] * fov_integral) / self.n_showers @@ -334,3 +329,69 @@ 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_mask(self, bin_edges_lon, bin_edges_lat, radius, subpixels=20): + # 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_mask = np.ones(mask_simple.shape) + area_mask[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_diameter_lon = bin_width_lon / np.sqrt(2) + bin_max_diameter_lat = bin_width_lat / np.sqrt(2) + + fov_edge_mask = np.logical_and( + radius_bin_center < radius + bin_max_diameter_lon, + radius_bin_center > radius - bin_max_diameter_lat, + ) + + # get indices of relevant bin corners + corner_idx = np.nonzero(fov_edge_mask) + + # define start and endpoints for subpixels + bin_grid_lon, bin_grid_lat = np.meshgrid(bin_edges_lon, bin_edges_lat) + 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 = radius_subpixel <= radius + + # calculates the fraction of subpixel centers within the FOV + FOV_covered_area = mask.sum(axis=(1,2)) / (subpixels ** 2) + + area_mask[corner_idx] = FOV_covered_area + + return area_mask \ No newline at end of file diff --git a/pyirf/tests/test_simulations.py b/pyirf/tests/test_simulations.py index fbb5ac335..649ec79f1 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 From 46746794dcc857cb0bd099ed9e4372266e9ce431 Mon Sep 17 00:00:00 2001 From: Luca Di Bella Date: Fri, 1 Mar 2024 16:48:47 +0100 Subject: [PATCH 33/48] Fix position angle function, add rectangle solid angle function and tests --- pyirf/tests/test_utils.py | 22 ++++++++++++++ pyirf/utils.py | 60 +++++++++++++++++++++++++++++++++++++-- 2 files changed, 79 insertions(+), 3 deletions(-) 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 4be243521..95fab4eb0 100644 --- a/pyirf/utils.py +++ b/pyirf/utils.py @@ -1,7 +1,7 @@ import numpy as np import astropy.units as u from astropy.coordinates import angular_separation -from astropy.coordinates import position_angle +from .coordinates import fov_coords_theta_phi, fov_coords_lon_lat from .compat import COPY_IF_NEEDED from .exceptions import MissingColumns, WrongColumnUnit @@ -109,14 +109,42 @@ def calculate_source_fov_position_angle(events, prefix="true"): Position angle of the true positions relative to the pointing positions in the sky. """ - phi = position_angle( + _, phi = 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 = fov_coords_lon_lat( events[f"{prefix}_az"], events[f"{prefix}_alt"], + events["pointing_az"], + events["pointing_alt"], ) - return phi.to(u.deg) + return lon.to(u.deg), lat.to(u.deg) def check_histograms(hist1, hist2, key="reco_energy"): @@ -159,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. From 5e3825410e18ccf79faf49c4c389bfa7e334d1ff Mon Sep 17 00:00:00 2001 From: Luca Di Bella Date: Fri, 1 Mar 2024 16:50:17 +0100 Subject: [PATCH 34/48] Update effective area 3D functions and add corresponding tests --- pyirf/irf/effective_area.py | 19 ++-- pyirf/irf/tests/test_effective_area.py | 145 +++++++++++++++++++++++++ 2 files changed, 157 insertions(+), 7 deletions(-) diff --git a/pyirf/irf/effective_area.py b/pyirf/irf/effective_area.py index 102a3add6..c076e33c2 100644 --- a/pyirf/irf/effective_area.py +++ b/pyirf/irf/effective_area.py @@ -140,7 +140,12 @@ def effective_area_3d_polar( def effective_area_3d_lonlat( - selected_events, simulation_info, energy_bins, fov_longitude_bins, fov_latitude_bins + selected_events, + simulation_info, + 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. @@ -164,22 +169,22 @@ def effective_area_3d_lonlat( area = np.pi * simulation_info.max_impact**2 hist_simulated = simulation_info.calculate_n_showers_3d_lonlat( - energy_bins, fov_longitude_bins, fov_latitude_bins + energy_bins, fov_longitude_bins, fov_latitude_bins, subpixels=subpixels ) - hist_selected, _ = np.histogramdd( - np.column_stack( + 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=( + ) + bins = ( 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 From 8e017c5a670683903d5ce34c41235b899d4267de Mon Sep 17 00:00:00 2001 From: Luca Di Bella Date: Mon, 15 Apr 2024 13:54:45 +0200 Subject: [PATCH 35/48] Fix code review issues --- pyirf/coordinates.py | 24 +++++++++++++++------ pyirf/irf/effective_area.py | 12 +++++------ pyirf/simulations.py | 38 ++++++++++++++++----------------- pyirf/tests/test_coordinates.py | 6 +++--- pyirf/tests/test_simulations.py | 2 +- pyirf/utils.py | 6 +++--- 6 files changed, 49 insertions(+), 39 deletions(-) diff --git a/pyirf/coordinates.py b/pyirf/coordinates.py index 96be28054..0c961d167 100644 --- a/pyirf/coordinates.py +++ b/pyirf/coordinates.py @@ -7,8 +7,12 @@ ] -def fov_coords_lon_lat(lon, lat, pointing_lon, pointing_lat): - """Transform sky coordinates to field-of-view longitude-latitude coordinates. +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 ---------- @@ -38,8 +42,12 @@ def fov_coords_lon_lat(lon, lat, pointing_lon, pointing_lat): return -target_fov.lon, target_fov.lat -def fov_coords_theta_phi(az, alt, pointing_az, pointing_alt): - """Transform sky coordinates to field-of-view theta-phi coordinates. +def gadf_fov_coords_theta_phi(az, alt, pointing_az, pointing_alt): + """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 ---------- @@ -51,12 +59,14 @@ def fov_coords_theta_phi(az, alt, pointing_az, pointing_alt): Returns ------- - lon_t, lat_t : `~astropy.units.Quantity` + theta, phi : `~astropy.units.Quantity` Transformed field-of-view coordinate. """ theta = angular_separation(pointing_az, pointing_alt, az, alt) - + + # astropy defines the position angle as increasing towards east of north phi = position_angle(pointing_az, pointing_alt, az, alt) - + + # 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/effective_area.py b/pyirf/irf/effective_area.py index c076e33c2..a204cd556 100644 --- a/pyirf/irf/effective_area.py +++ b/pyirf/irf/effective_area.py @@ -92,7 +92,7 @@ def effective_area_per_energy_and_fov( def effective_area_3d_polar( selected_events, simulation_info, - energy_bins, + true_energy_bins, fov_offset_bins, fov_position_angle_bins, ): @@ -118,7 +118,7 @@ def effective_area_3d_polar( area = np.pi * simulation_info.max_impact**2 hist_simulated = simulation_info.calculate_n_showers_3d_polar( - energy_bins, fov_offset_bins, fov_position_angle_bins + true_energy_bins, fov_offset_bins, fov_position_angle_bins ) hist_selected, _ = np.histogramdd( @@ -130,7 +130,7 @@ def effective_area_3d_polar( ] ), bins=( - energy_bins.to_value(u.TeV), + true_energy_bins.to_value(u.TeV), fov_offset_bins.to_value(u.deg), fov_position_angle_bins.to_value(u.rad), ), @@ -142,7 +142,7 @@ def effective_area_3d_polar( def effective_area_3d_lonlat( selected_events, simulation_info, - energy_bins, + true_energy_bins, fov_longitude_bins, fov_latitude_bins, subpixels=20, @@ -169,7 +169,7 @@ def effective_area_3d_lonlat( area = np.pi * simulation_info.max_impact**2 hist_simulated = simulation_info.calculate_n_showers_3d_lonlat( - energy_bins, fov_longitude_bins, fov_latitude_bins, subpixels=subpixels + true_energy_bins, fov_longitude_bins, fov_latitude_bins, subpixels=subpixels ) selected_columns = np.column_stack( @@ -180,7 +180,7 @@ def effective_area_3d_lonlat( ] ) bins = ( - energy_bins.to_value(u.TeV), + true_energy_bins.to_value(u.TeV), fov_longitude_bins.to_value(u.deg), fov_latitude_bins.to_value(u.deg), ) diff --git a/pyirf/simulations.py b/pyirf/simulations.py index fc99bdc78..63defbd71 100644 --- a/pyirf/simulations.py +++ b/pyirf/simulations.py @@ -210,10 +210,10 @@ def calculate_n_showers_3d_polar( 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 @@ -250,8 +250,8 @@ def calculate_n_showers_3d_lonlat( The actual numbers will follow a poissionian distribution around this expected value. """ - - fov_mask = _fov_lonlat_grid_overlap_mask( + + fov_overlap = _fov_lonlat_grid_overlap( self, fov_longitude_bins, fov_latitude_bins, self.viewcone_max, subpixels=subpixels ) @@ -269,7 +269,7 @@ def calculate_n_showers_3d_lonlat( fov_integral = shower_density * bin_area e_integral = self.calculate_n_showers_per_energy(energy_bins) - fov_integral = fov_mask * fov_integral + fov_integral = fov_overlap * fov_integral return (e_integral[:, np.newaxis, np.newaxis] * fov_integral) / self.n_showers @@ -330,7 +330,7 @@ def _viewcone_pdf_integral(viewcone_min, viewcone_max, fov_low, fov_high): return np.squeeze(integral) return integral -def _fov_lonlat_grid_overlap_mask(self, bin_edges_lon, bin_edges_lat, radius, subpixels=20): +def _fov_lonlat_grid_overlap(self, bin_edges_lon, bin_edges_lat, radius, subpixels=20): # define grid of bin centers fov_bin_centers_lon = bin_center(bin_edges_lon) fov_bin_centers_lat = bin_center(bin_edges_lat) @@ -338,7 +338,7 @@ def _fov_lonlat_grid_overlap_mask(self, bin_edges_lon, bin_edges_lat, radius, su 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) @@ -346,23 +346,23 @@ def _fov_lonlat_grid_overlap_mask(self, bin_edges_lon, bin_edges_lat, radius, su mask_simple = np.logical_or( radius_bin_center > self.viewcone_max, radius_bin_center < self.viewcone_min ) - area_mask = np.ones(mask_simple.shape) - area_mask[mask_simple] = 0 + 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_diameter_lon = bin_width_lon / np.sqrt(2) bin_max_diameter_lat = bin_width_lat / np.sqrt(2) - + fov_edge_mask = np.logical_and( radius_bin_center < radius + bin_max_diameter_lon, radius_bin_center > radius - bin_max_diameter_lat, ) - - # get indices of relevant bin corners + + # get indices of relevant bin corners corner_idx = np.nonzero(fov_edge_mask) - + # define start and endpoints for subpixels bin_grid_lon, bin_grid_lat = np.meshgrid(bin_edges_lon, bin_edges_lat) edges_lon = np.array( @@ -371,9 +371,9 @@ def _fov_lonlat_grid_overlap_mask(self, bin_edges_lon, bin_edges_lat, radius, su 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) @@ -384,7 +384,7 @@ def _fov_lonlat_grid_overlap_mask(self, bin_edges_lon, bin_edges_lat, radius, su ) 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 = radius_subpixel <= radius @@ -392,6 +392,6 @@ def _fov_lonlat_grid_overlap_mask(self, bin_edges_lon, bin_edges_lat, radius, su # calculates the fraction of subpixel centers within the FOV FOV_covered_area = mask.sum(axis=(1,2)) / (subpixels ** 2) - area_mask[corner_idx] = FOV_covered_area - - return area_mask \ No newline at end of file + 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 index 5a04a5d19..5196d6d7b 100644 --- a/pyirf/tests/test_coordinates.py +++ b/pyirf/tests/test_coordinates.py @@ -40,19 +40,19 @@ def test_fov_coords_theta_phi(): ) assert u.isclose(theta, 1 * u.deg) assert u.isclose(phi, 0 * u.deg) - + theta, phi = fov_coords_theta_phi( alt=-1 * u.deg, az=0 * u.deg, pointing_alt=0 * u.deg, pointing_az=0 * u.deg ) assert u.isclose(theta, 1 * u.deg) assert u.isclose(phi, 180 * u.deg) - + theta, phi = fov_coords_theta_phi( alt=0 * u.deg, az=-1 * u.deg, pointing_alt=0 * u.deg, pointing_az=0 * u.deg ) assert u.isclose(theta, 1 * u.deg) assert u.isclose(phi, 90 * u.deg) - + theta, phi = fov_coords_theta_phi( alt=0 * u.deg, az=1 * u.deg, pointing_alt=0 * u.deg, pointing_az=0 * u.deg ) diff --git a/pyirf/tests/test_simulations.py b/pyirf/tests/test_simulations.py index 649ec79f1..5ccabf827 100644 --- a/pyirf/tests/test_simulations.py +++ b/pyirf/tests/test_simulations.py @@ -169,7 +169,7 @@ def test_integrate_3d_lonlat(): 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) diff --git a/pyirf/utils.py b/pyirf/utils.py index 95fab4eb0..39ddd83f3 100644 --- a/pyirf/utils.py +++ b/pyirf/utils.py @@ -189,7 +189,7 @@ def cone_solid_angle(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] @@ -200,7 +200,7 @@ def rectangle_solid_angle(lon_low, lon_high, lat_low, lat_high): 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] @@ -208,7 +208,7 @@ def rectangle_solid_angle(lon_low, lon_high, lat_low, lat_high): """ 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 From 2a91a87615ff490b5f4c08cf4f64b2bfec830420 Mon Sep 17 00:00:00 2001 From: Luca Di Bella Date: Thu, 18 Apr 2024 10:47:33 +0200 Subject: [PATCH 36/48] Update function names in __all__ --- pyirf/coordinates.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pyirf/coordinates.py b/pyirf/coordinates.py index 0c961d167..425850bcf 100644 --- a/pyirf/coordinates.py +++ b/pyirf/coordinates.py @@ -2,8 +2,8 @@ from astropy.coordinates import SkyCoord, SkyOffsetFrame, angular_separation, position_angle __all__ = [ - 'fov_coords_lon_lat', - 'fov_coords_theta_phi' + 'gadf_fov_coords_lon_lat', + 'gadf_fov_coords_theta_phi', ] From b516419104e2a717516de008eb8d59d70e6d9a58 Mon Sep 17 00:00:00 2001 From: Luca Di Bella Date: Thu, 18 Apr 2024 11:04:19 +0200 Subject: [PATCH 37/48] Fix all changed occurences of gadf coordinate function calls --- pyirf/tests/test_coordinates.py | 24 ++++++++++++------------ pyirf/utils.py | 6 +++--- 2 files changed, 15 insertions(+), 15 deletions(-) diff --git a/pyirf/tests/test_coordinates.py b/pyirf/tests/test_coordinates.py index 5196d6d7b..b687a9d53 100644 --- a/pyirf/tests/test_coordinates.py +++ b/pyirf/tests/test_coordinates.py @@ -1,18 +1,18 @@ from numpy.testing import assert_allclose import astropy.units as u -def test_fov_coords_lon_lat(): - from pyirf.coordinates import fov_coords_lon_lat +def test_gadf_fov_coords_lon_lat(): + from pyirf.coordinates import gadf_fov_coords_lon_lat # test some simple cases - lon, lat = fov_coords_lon_lat(1 * u.deg, 1 * u.deg, 0 * u.deg, 0 * u.deg) + 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 = fov_coords_lon_lat(269 * u.deg, 0 * u.deg, 270 * u.deg, 0 * u.deg) + 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 = fov_coords_lon_lat(1 * u.deg, 60 * u.deg, 0 * u.deg, 60 * u.deg) + 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) @@ -22,7 +22,7 @@ def test_fov_coords_lon_lat(): 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 = fov_coords_lon_lat( + lon, lat = gadf_fov_coords_lon_lat( az * u.deg, alt * u.deg, az_pointing * u.deg, alt_pointing * u.deg ) assert_allclose( @@ -32,28 +32,28 @@ def test_fov_coords_lon_lat(): lat.value, [-1.60829115, -1.19643974, 0.45800984, 3.26844192], rtol=1e-5 ) -def test_fov_coords_theta_phi(): - from pyirf.coordinates import fov_coords_theta_phi +def test_gadf_fov_coords_theta_phi(): + from pyirf.coordinates import gadf_fov_coords_theta_phi - theta, phi = fov_coords_theta_phi( + theta, phi = gadf_fov_coords_theta_phi( alt=1 * u.deg, az=0 * u.deg, pointing_alt=0 * u.deg, pointing_az=0 * u.deg ) assert u.isclose(theta, 1 * u.deg) assert u.isclose(phi, 0 * u.deg) - theta, phi = fov_coords_theta_phi( + theta, phi = gadf_fov_coords_theta_phi( alt=-1 * u.deg, az=0 * u.deg, pointing_alt=0 * u.deg, pointing_az=0 * u.deg ) assert u.isclose(theta, 1 * u.deg) assert u.isclose(phi, 180 * u.deg) - theta, phi = fov_coords_theta_phi( + theta, phi = gadf_fov_coords_theta_phi( alt=0 * u.deg, az=-1 * u.deg, pointing_alt=0 * u.deg, pointing_az=0 * u.deg ) assert u.isclose(theta, 1 * u.deg) assert u.isclose(phi, 90 * u.deg) - theta, phi = fov_coords_theta_phi( + theta, phi = gadf_fov_coords_theta_phi( alt=0 * u.deg, az=1 * u.deg, pointing_alt=0 * u.deg, pointing_az=0 * u.deg ) assert u.isclose(theta, 1 * u.deg) diff --git a/pyirf/utils.py b/pyirf/utils.py index 39ddd83f3..e1d41c631 100644 --- a/pyirf/utils.py +++ b/pyirf/utils.py @@ -1,7 +1,7 @@ import numpy as np import astropy.units as u from astropy.coordinates import angular_separation -from .coordinates import fov_coords_theta_phi, fov_coords_lon_lat +from .coordinates import gadf_fov_coords_theta_phi, gadf_fov_coords_lon_lat from .compat import COPY_IF_NEEDED from .exceptions import MissingColumns, WrongColumnUnit @@ -109,7 +109,7 @@ def calculate_source_fov_position_angle(events, prefix="true"): Position angle of the true positions relative to the pointing positions in the sky. """ - _, phi = fov_coords_theta_phi( + _, phi = gadf_fov_coords_theta_phi( events[f"{prefix}_az"], events[f"{prefix}_alt"], events["pointing_az"], @@ -137,7 +137,7 @@ def calculate_source_fov_lonlat(events, prefix="true"): Position angle of the true positions relative to the pointing positions in the sky. """ - lon, lat = fov_coords_lon_lat( + lon, lat = gadf_fov_coords_lon_lat( events[f"{prefix}_az"], events[f"{prefix}_alt"], events["pointing_az"], From e18e1aeb138a93e96e62130cc447882731fe8eb3 Mon Sep 17 00:00:00 2001 From: Luca Di Bella Date: Tue, 30 Apr 2024 15:26:06 +0200 Subject: [PATCH 38/48] Change az/alt to lon/lat due to support of both AltAz and ICRS --- pyirf/coordinates.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/pyirf/coordinates.py b/pyirf/coordinates.py index 425850bcf..d8e2eca66 100644 --- a/pyirf/coordinates.py +++ b/pyirf/coordinates.py @@ -42,7 +42,7 @@ def gadf_fov_coords_lon_lat(lon, lat, pointing_lon, pointing_lat): return -target_fov.lon, target_fov.lat -def gadf_fov_coords_theta_phi(az, alt, pointing_az, pointing_alt): +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. @@ -51,9 +51,9 @@ def gadf_fov_coords_theta_phi(az, alt, pointing_az, pointing_alt): Parameters ---------- - az, alt : `~astropy.units.Quantity` + lon, lat : `~astropy.units.Quantity` Sky coordinate to be transformed. - pointing_az, pointing_alt : `~astropy.units.Quantity` + pointing_lon, pointing_lat : `~astropy.units.Quantity` Coordinate specifying the pointing position. (i.e. the center of the field of view.) @@ -63,10 +63,10 @@ def gadf_fov_coords_theta_phi(az, alt, pointing_az, pointing_alt): Transformed field-of-view coordinate. """ - theta = angular_separation(pointing_az, pointing_alt, az, alt) + theta = angular_separation(pointing_lon, pointing_lat, lon, lat) # astropy defines the position angle as increasing towards east of north - phi = position_angle(pointing_az, pointing_alt, az, alt) + 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 From 3162c4a235ff7bb83c4bb9b8bc4f0baa72b7ae45 Mon Sep 17 00:00:00 2001 From: Luca Di Bella Date: Tue, 30 Apr 2024 15:43:32 +0200 Subject: [PATCH 39/48] Add newsfragment --- docs/changes/281.feature.rst | 2 ++ 1 file changed, 2 insertions(+) create mode 100644 docs/changes/281.feature.rst 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. From 8cd0e756170c7dd40e51841cb9552ce447d3b228 Mon Sep 17 00:00:00 2001 From: Luca Di Bella Date: Thu, 2 May 2024 16:45:40 +0200 Subject: [PATCH 40/48] Properly handle viewcone_min > 0, add docstring --- pyirf/simulations.py | 66 +++++++++++++++++++++++++++------ pyirf/tests/test_coordinates.py | 8 ++-- 2 files changed, 59 insertions(+), 15 deletions(-) diff --git a/pyirf/simulations.py b/pyirf/simulations.py index 63defbd71..d03346d4d 100644 --- a/pyirf/simulations.py +++ b/pyirf/simulations.py @@ -252,7 +252,7 @@ def calculate_n_showers_3d_lonlat( """ fov_overlap = _fov_lonlat_grid_overlap( - self, fov_longitude_bins, fov_latitude_bins, self.viewcone_max, subpixels=subpixels + self, fov_longitude_bins, fov_latitude_bins, subpixels=subpixels ) bin_grid_lon, bin_grid_lat = np.meshgrid(fov_longitude_bins,fov_latitude_bins) @@ -330,7 +330,37 @@ def _viewcone_pdf_integral(viewcone_min, viewcone_max, fov_low, fov_high): return np.squeeze(integral) return integral -def _fov_lonlat_grid_overlap(self, bin_edges_lon, bin_edges_lat, radius, subpixels=20): +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.all(bin_dist <= self.viewcone_min.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) @@ -340,31 +370,42 @@ def _fov_lonlat_grid_overlap(self, bin_edges_lon, bin_edges_lat, radius, subpixe ) # calculate angular separation of bin centers to FOV center - radius_bin_center = angular_separation(bin_centers_grid_lon, bin_centers_grid_lat, 0, 0) + 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 + 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_diameter_lon = bin_width_lon / np.sqrt(2) - bin_max_diameter_lat = bin_width_lat / np.sqrt(2) + bin_max_radius = np.sqrt(bin_width_lon ** 2 + bin_width_lat ** 2) * 0.5 - fov_edge_mask = np.logical_and( - radius_bin_center < radius + bin_max_diameter_lon, - radius_bin_center > radius - bin_max_diameter_lat, + 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 - bin_grid_lon, bin_grid_lat = np.meshgrid(bin_edges_lon, bin_edges_lat) edges_lon = np.array( [bin_grid_lon[corner_idx], bin_grid_lon[corner_idx] + bin_width_lon] ) @@ -387,7 +428,10 @@ def _fov_lonlat_grid_overlap(self, bin_edges_lon, bin_edges_lat, radius, subpixe # make mask with subpixels inside the FOV radius_subpixel = angular_separation(subpixel_grid_lon, subpixel_grid_lat, 0, 0) - mask = radius_subpixel <= radius + 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) diff --git a/pyirf/tests/test_coordinates.py b/pyirf/tests/test_coordinates.py index b687a9d53..e96a639dd 100644 --- a/pyirf/tests/test_coordinates.py +++ b/pyirf/tests/test_coordinates.py @@ -36,25 +36,25 @@ def test_gadf_fov_coords_theta_phi(): from pyirf.coordinates import gadf_fov_coords_theta_phi theta, phi = gadf_fov_coords_theta_phi( - alt=1 * u.deg, az=0 * u.deg, pointing_alt=0 * u.deg, pointing_az=0 * u.deg + 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( - alt=-1 * u.deg, az=0 * u.deg, pointing_alt=0 * u.deg, pointing_az=0 * u.deg + 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( - alt=0 * u.deg, az=-1 * u.deg, pointing_alt=0 * u.deg, pointing_az=0 * u.deg + 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( - alt=0 * u.deg, az=1 * u.deg, pointing_alt=0 * u.deg, pointing_az=0 * u.deg + 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 From 82b60b37836c82cf14c8de38e8cb56b8d6caa412 Mon Sep 17 00:00:00 2001 From: Luca Di Bella Date: Tue, 14 May 2024 12:06:35 +0200 Subject: [PATCH 41/48] Add missing condition to all_outside check --- pyirf/simulations.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/pyirf/simulations.py b/pyirf/simulations.py index d03346d4d..d7f6c3be9 100644 --- a/pyirf/simulations.py +++ b/pyirf/simulations.py @@ -349,7 +349,10 @@ def _fov_lonlat_grid_overlap(self, bin_edges_lon, bin_edges_lat, subpixels=20): 0, ).to_value(u.deg) - all_outside = np.all(bin_dist <= self.viewcone_min.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)), From a63a56bca89b664302de55c4aaf330d86032b15f Mon Sep 17 00:00:00 2001 From: Michele Peresano Date: Tue, 12 Dec 2023 11:58:52 +0100 Subject: [PATCH 42/48] Update energy_dispersion_to_migration Remove stale debugging code Add news fragment Update energy_dispersion_to_migration - divide by the migration bin widths before resampling - minor style and format changes Fix pdf conversion to probability, remove final norm --- docs/changes/273.bugfix.rst | 4 ++++ pyirf/irf/energy_dispersion.py | 34 +++++++++++++++++++++------------- 2 files changed, 25 insertions(+), 13 deletions(-) create mode 100644 docs/changes/273.bugfix.rst diff --git a/docs/changes/273.bugfix.rst b/docs/changes/273.bugfix.rst new file mode 100644 index 000000000..607ffa71f --- /dev/null +++ b/docs/changes/273.bugfix.rst @@ -0,0 +1,4 @@ +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. \ No newline at end of file diff --git a/pyirf/irf/energy_dispersion.py b/pyirf/irf/energy_dispersion.py index de631e426..f6100adda 100644 --- a/pyirf/irf/energy_dispersion.py +++ b/pyirf/irf/energy_dispersion.py @@ -29,7 +29,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. @@ -88,8 +91,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 +113,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 +144,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( From f3b4d04da5607f7415cdc46b168e0829821f1edb Mon Sep 17 00:00:00 2001 From: Michele Peresano Date: Tue, 12 Dec 2023 14:53:06 +0100 Subject: [PATCH 43/48] Add function to compute energy migration matrix from events Simplify normalization code for migra matrix --- pyirf/irf/energy_dispersion.py | 53 +++++++++++++++++++++++++++++++++- 1 file changed, 52 insertions(+), 1 deletion(-) diff --git a/pyirf/irf/energy_dispersion.py b/pyirf/irf/energy_dispersion.py index f6100adda..c3d7bcac8 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", ] @@ -83,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 true + energy acis 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, From 7f2d45d2bbfeefb0da498cd3e6886bfe08c40dcd Mon Sep 17 00:00:00 2001 From: Michele Peresano Date: Tue, 12 Dec 2023 15:03:55 +0100 Subject: [PATCH 44/48] Add unit test for energy migration matrix from events update unit test --- pyirf/irf/tests/test_energy_dispersion.py | 104 +++++++++++++++++----- 1 file changed, 80 insertions(+), 24 deletions(-) diff --git a/pyirf/irf/tests/test_energy_dispersion.py b/pyirf/irf/tests/test_energy_dispersion.py index 036d106de..ec0e444fb 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 @@ -133,14 +143,57 @@ def test_energy_dispersion_to_migration(): assert migration_matrix.shape[2] == dispersion_matrix.shape[2] # test that all migrations are included for central energies - assert np.isclose(migration_matrix.sum(axis=1).max(), 1, rtol=0.01) + np.allclose(migration_matrix.sum(axis=1).max(), 1, rtol=0.1) # 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( From 2bd128c57990dda57b087c43928449069c57fecd Mon Sep 17 00:00:00 2001 From: Michele Peresano Date: Tue, 12 Dec 2023 15:10:23 +0100 Subject: [PATCH 45/48] Update test_energy_dispersion_to_migration Fix test_energy_dispersion_to_migration --- pyirf/irf/tests/test_energy_dispersion.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyirf/irf/tests/test_energy_dispersion.py b/pyirf/irf/tests/test_energy_dispersion.py index ec0e444fb..451511ffc 100644 --- a/pyirf/irf/tests/test_energy_dispersion.py +++ b/pyirf/irf/tests/test_energy_dispersion.py @@ -143,7 +143,7 @@ def test_energy_dispersion_to_migration(): assert migration_matrix.shape[2] == dispersion_matrix.shape[2] # test that all migrations are included for central energies - np.allclose(migration_matrix.sum(axis=1).max(), 1, rtol=0.1) + assert np.isclose(migration_matrix.sum(axis=1).max(), 1, rtol=0.01) # test that migrations dont always sum to 1 (since some energies are # not included in the matrix) From 4bc8f5780b13e4d13f0c72a13e793d590f7ebf11 Mon Sep 17 00:00:00 2001 From: Michele Peresano Date: Thu, 7 Nov 2024 11:28:39 +0100 Subject: [PATCH 46/48] fix: docstring irf.energy_dispersion.energy_migration_matrix --- pyirf/irf/energy_dispersion.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pyirf/irf/energy_dispersion.py b/pyirf/irf/energy_dispersion.py index c3d7bcac8..f34c3546d 100644 --- a/pyirf/irf/energy_dispersion.py +++ b/pyirf/irf/energy_dispersion.py @@ -103,8 +103,8 @@ def energy_migration_matrix( Returns ------- matrix : array-like - Migration matrix as probabilities along the true - energy acis with shape + 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. """ From 82e20354fcba68b207ec5d223c7ada4e2bd65a8f Mon Sep 17 00:00:00 2001 From: Tomas Bylund Date: Tue, 16 Jan 2024 14:50:52 +0100 Subject: [PATCH 47/48] Add basic facilities for 3D background irfs --- pyirf/io/gadf.py | 51 ++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 51 insertions(+) diff --git a/pyirf/io/gadf.py b/pyirf/io/gadf.py index 963b6de78..335d5146c 100644 --- a/pyirf/io/gadf.py +++ b/pyirf/io/gadf.py @@ -318,6 +318,57 @@ def create_background_3d_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_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, From 03b96264f9e2c47863d0c6ed0abfa9519697a0e0 Mon Sep 17 00:00:00 2001 From: Tomas Bylund Date: Fri, 8 Nov 2024 19:27:14 +0100 Subject: [PATCH 48/48] Fixes --- docs/changes/276.feature.rst | 1 + pyirf/io/gadf.py | 11 +++++--- pyirf/irf/__init__.py | 4 +-- pyirf/irf/background.py | 40 ++++++++++++++---------------- pyirf/irf/tests/test_background.py | 13 +++++----- 5 files changed, 36 insertions(+), 33 deletions(-) create mode 100644 docs/changes/276.feature.rst 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/pyirf/io/gadf.py b/pyirf/io/gadf.py index 335d5146c..764fc8118 100644 --- a/pyirf/io/gadf.py +++ b/pyirf/io/gadf.py @@ -268,7 +268,7 @@ def create_background_3d_hdu( reco_energy_bins, fov_lon_bins, fov_lat_bins, - alignment = "ALTAZ", + alignment="ALTAZ", extname="BACKGROUND", **header_cards, ): @@ -289,7 +289,7 @@ def create_background_3d_hdu( fov_lat_bins: astropy.units.Quantity[angle] Bin edges in the field of view system, becomes the DETY values alignment: str - Wheter the FOV coordinates are aligned with the ALTAZ or RADEC system, more details at + 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 @@ -310,8 +310,11 @@ def create_background_3d_hdu( header["HDUCLAS1"] = "RESPONSE" header["HDUCLAS2"] = "BKG" header["HDUCLAS3"] = "FULL-ENCLOSURE" - header["HDUCLAS4"] = "BKG_2D" - header["FOVALIGN"] = alignment + 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) diff --git a/pyirf/irf/__init__.py b/pyirf/irf/__init__.py index b002b4289..82bae2372 100644 --- a/pyirf/irf/__init__.py +++ b/pyirf/irf/__init__.py @@ -7,7 +7,7 @@ ) from .energy_dispersion import energy_dispersion from .psf import psf_table -from .background import background_2d, background_3d +from .background import background_2d, background_3d_lonlat __all__ = [ "effective_area", @@ -18,5 +18,5 @@ "energy_dispersion", "psf_table", "background_2d", - "background_3d", + "background_3d_lonlat", ] diff --git a/pyirf/irf/background.py b/pyirf/irf/background.py index a5f9b3b0d..55b41f0ee 100644 --- a/pyirf/irf/background.py +++ b/pyirf/irf/background.py @@ -58,7 +58,7 @@ def background_2d(events, reco_energy_bins, fov_offset_bins, t_obs): return bg_rate.to(BACKGROUND_UNIT) -def background_3d(events, reco_energy_bins, fov_offset_bins, t_obs): +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. @@ -73,8 +73,12 @@ def background_3d(events, reco_energy_bins, fov_offset_bins, t_obs): `reco_energy`, `weight`. reco_energy: astropy.units.Quantity[energy] The bins in reconstructed energy to be used for the IRF - fov_offset_bins: astropy.units.Quantity[angle] - The bins in the field of view offset to be used for the IRF, either a (N,) or (1,N) or a (2,N) array + 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. @@ -85,24 +89,16 @@ def background_3d(events, reco_energy_bins, fov_offset_bins, t_obs): The background rate as particles per energy, time and solid angle in the specified bins. - Shape: (len(reco_energy_bins) - 1, len(fov_offset_bins) - 1, len(fov_offset_bins) - 1) + Shape: (len(reco_energy_bins) - 1, len(fov_lon_bins) - 1, len(fov_lon_bins) - 1) """ - if (fov_offset_bins.shape[0] == 1) or (len(fov_offset_bins.shape) == 1): - fov_x_offset_bins = fov_offset_bins - fov_y_offset_bins = fov_offset_bins - elif fov_offset_bins.shape[0] == 2: - fov_x_offset_bins = fov_offset_bins[0, :] - fov_y_offset_bins = fov_offset_bins[1, :] - else: - raise ValueError( - f"fov_offset_bins must be eiher (N,) or (1,N) or (2,N), found {fov_offset_bins.shape}" - ) + 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), + events["reco_fov_lon"].to_value(u.deg), + events["reco_fov_lat"].to_value(u.deg), ], bins=[ reco_energy_bins.to_value(u.TeV), @@ -115,12 +111,14 @@ def background_3d(events, reco_energy_bins, fov_offset_bins, t_obs): # 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.T / bin_width_energy).T + 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:]) + 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/tests/test_background.py b/pyirf/irf/tests/test_background.py index 5d3e009ae..cfb4e2b18 100644 --- a/pyirf/irf/tests/test_background.py +++ b/pyirf/irf/tests/test_background.py @@ -41,8 +41,8 @@ def test_background(): ) -def test_background_3d(): - from pyirf.irf import background_3d +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 @@ -103,10 +103,11 @@ def test_background_3d(): } ) - bkg_rate = background_3d( + bkg_rate = background_3d_lonlat( selected_events, reco_energy_bins=reco_energy_bins, - fov_offset_bins=np.vstack((fov_lat_bins, fov_lon_bins)), + fov_lon_bins=fov_lon_bins, + fov_lat_bins=fov_lat_bins, t_obs=t_obs, ) assert bkg_rate.shape == ( @@ -118,11 +119,11 @@ def test_background_3d(): # Convert to counts, project to energy axis, and check counts round-trip correctly assert np.allclose( - (bin_solid_angle * (bkg_rate.T * bin_width_energy).T).sum(axis=(1, 2)) * t_obs, + (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.T * bin_width_energy).T).sum(axis=(0, 1)) * t_obs, + (bin_solid_angle * bkg_rate * bin_width_energy[:, np.newaxis, np.newaxis]).sum(axis=(0, 1)) * t_obs, 2 * [N_tot // 2], )