From ab98804dbf79c89e35aea3b62b3982967bf0305f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Bruno=20Kh=C3=A9lifi?= Date: Mon, 14 Oct 2024 14:34:09 +0200 Subject: [PATCH 01/55] new --- docs/_templates/custom-footer.html | 4 ++++ 1 file changed, 4 insertions(+) create mode 100644 docs/_templates/custom-footer.html diff --git a/docs/_templates/custom-footer.html b/docs/_templates/custom-footer.html new file mode 100644 index 0000000000..84facf5f1b --- /dev/null +++ b/docs/_templates/custom-footer.html @@ -0,0 +1,4 @@ +{% block extrafooter %} + {{super}} + Data Privacy. +{% endblock %} From 38cba1603cfc3d3c7c8a4b6faedea012d24502cc Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Bruno=20Kh=C3=A9lifi?= Date: Mon, 14 Oct 2024 14:34:44 +0200 Subject: [PATCH 02/55] update footer --- docs/conf.py | 15 +++++++-------- 1 file changed, 7 insertions(+), 8 deletions(-) diff --git a/docs/conf.py b/docs/conf.py index 1a91a01516..5d2ed8d6f8 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -112,7 +112,7 @@ def setup(app): # List of patterns, relative to source directory, that match files and # directories to ignore when looking for source files. -exclude_patterns.append("_templates") +#exclude_patterns.append("_templates") exclude_patterns.append("**.ipynb_checkpoints") exclude_patterns.append("user-guide/model-gallery/*/*.ipynb") exclude_patterns.append("user-guide/model-gallery/*/*.md5") @@ -182,6 +182,11 @@ def setup(app): # Static files to copy after template files html_static_path = ["_static"] +html_css_files = ["custom.css"] +html_js_files = ["matomo.js"] + +templates_path = ["_templates"] + # The name of an image file (within the static path) to use as favicon of the # docs. This file should be a Windows icon file (.ico) being 16x16 or 32x32 @@ -236,7 +241,7 @@ def setup(app): "navbar_end": ["version-switcher", "theme-switcher", "navbar-icon-links"], "navigation_with_keys": True, # footers - "footer_start": ["copyright"], + "footer_start": ["copyright","custom-footer.html"], "footer_center": ["last-updated"], "footer_end": ["sphinx-version", "theme-version"] } @@ -332,13 +337,7 @@ def setup(app): }, } -html_static_path = ["_static"] -html_css_files = ["custom.css"] -html_js_files = ["matomo.js"] - html_context = { "default_mode": "light", } -# Add-on to insert the Matomo tracker -templates_path = ['_templates'] From a49d7db71ab1691614eddc731f65e2ad7f635a35 Mon Sep 17 00:00:00 2001 From: Simone Date: Tue, 18 Jul 2023 11:55:54 +0200 Subject: [PATCH 03/55] use sky_to_fov also for RADEC FoVAlignment --- gammapy/makers/utils.py | 27 +++++++++++++++++++++++++++ 1 file changed, 27 insertions(+) diff --git a/gammapy/makers/utils.py b/gammapy/makers/utils.py index a7fd48e80f..cb2309b193 100644 --- a/gammapy/makers/utils.py +++ b/gammapy/makers/utils.py @@ -258,6 +258,33 @@ def make_map_background_irf( ) coords["energy"] = broadcast_axis_values_to_geom(geom, "energy", False) + # for backwards compatibility, obstime should be required + if obstime is None: + warnings.warn( + "Future versions of gammapy will require the obstime keyword for this function", + DeprecationWarning, + ) + obstime = pointing.obstime + + pointing_altaz = pointing.get_altaz(obstime) + altaz_coord = sky_coord.transform_to(pointing_altaz.frame) + + # Compute FOV coordinates of map relative to pointing + fov_lon, fov_lat = sky_to_fov( + altaz_coord.az, altaz_coord.alt, pointing_altaz.az, pointing_altaz.alt + ) + elif bkg.fov_alignment == FoVAlignment.RADEC: + fov_lon, fov_lat = sky_to_fov( + sky_coord.ra, sky_coord.dec, pointing_icrs.ra, pointing_icrs.dec + ) + else: + raise ValueError( + f"Unsupported background coordinate system: {bkg.fov_alignment!r}" + ) + + coords["fov_lon"] = fov_lon + coords["fov_lat"] = fov_lat + bkg_de = bkg.integrate_log_log(**coords, axis_name="energy") data = (bkg_de * d_omega * ontime).to_value("") From 3d840d5cbcbf9af0eeddc721019b0a61d90eb667 Mon Sep 17 00:00:00 2001 From: Quentin Remy Date: Wed, 17 Jul 2024 17:01:14 +0200 Subject: [PATCH 04/55] fix after rebase Signed-off-by: Quentin Remy --- gammapy/makers/utils.py | 37 ++++--------------------------------- 1 file changed, 4 insertions(+), 33 deletions(-) diff --git a/gammapy/makers/utils.py b/gammapy/makers/utils.py index cb2309b193..9b688dec29 100644 --- a/gammapy/makers/utils.py +++ b/gammapy/makers/utils.py @@ -3,7 +3,7 @@ import warnings import numpy as np import astropy.units as u -from astropy.coordinates import Angle, SkyOffsetFrame +from astropy.coordinates import Angle from astropy.table import Table from gammapy.data import FixedPointingInfo from gammapy.irf import BackgroundIRF, EDispMap, FoVAlignment, PSFMap @@ -83,11 +83,9 @@ def _get_fov_coords(pointing, irf, geom, use_region_center=True, obstime=None): altaz_coord.az, altaz_coord.alt, pointing_altaz.az, pointing_altaz.alt ) elif irf.fov_alignment == FoVAlignment.RADEC: - # Create OffsetFrame - frame = SkyOffsetFrame(origin=pointing_icrs) - pseudo_fov_coord = sky_coord.transform_to(frame) - fov_lon = pseudo_fov_coord.lon - fov_lat = pseudo_fov_coord.lat + fov_lon, fov_lat = sky_to_fov( + sky_coord.ra, sky_coord.dec, pointing_icrs.ra, pointing_icrs.dec + ) else: raise ValueError( f"Unsupported background coordinate system: {irf.fov_alignment!r}" @@ -258,33 +256,6 @@ def make_map_background_irf( ) coords["energy"] = broadcast_axis_values_to_geom(geom, "energy", False) - # for backwards compatibility, obstime should be required - if obstime is None: - warnings.warn( - "Future versions of gammapy will require the obstime keyword for this function", - DeprecationWarning, - ) - obstime = pointing.obstime - - pointing_altaz = pointing.get_altaz(obstime) - altaz_coord = sky_coord.transform_to(pointing_altaz.frame) - - # Compute FOV coordinates of map relative to pointing - fov_lon, fov_lat = sky_to_fov( - altaz_coord.az, altaz_coord.alt, pointing_altaz.az, pointing_altaz.alt - ) - elif bkg.fov_alignment == FoVAlignment.RADEC: - fov_lon, fov_lat = sky_to_fov( - sky_coord.ra, sky_coord.dec, pointing_icrs.ra, pointing_icrs.dec - ) - else: - raise ValueError( - f"Unsupported background coordinate system: {bkg.fov_alignment!r}" - ) - - coords["fov_lon"] = fov_lon - coords["fov_lat"] = fov_lat - bkg_de = bkg.integrate_log_log(**coords, axis_name="energy") data = (bkg_de * d_omega * ontime).to_value("") From b4b6a683cad6b6072b7d5ebb684a3f4b6d021001 Mon Sep 17 00:00:00 2001 From: Quentin Remy Date: Wed, 17 Jul 2024 17:10:31 +0200 Subject: [PATCH 05/55] icrs conversion Signed-off-by: Quentin Remy --- gammapy/makers/utils.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/gammapy/makers/utils.py b/gammapy/makers/utils.py index 9b688dec29..17b3476a43 100644 --- a/gammapy/makers/utils.py +++ b/gammapy/makers/utils.py @@ -84,7 +84,10 @@ def _get_fov_coords(pointing, irf, geom, use_region_center=True, obstime=None): ) elif irf.fov_alignment == FoVAlignment.RADEC: fov_lon, fov_lat = sky_to_fov( - sky_coord.ra, sky_coord.dec, pointing_icrs.ra, pointing_icrs.dec + sky_coord.icrs.ra, + sky_coord.icrs.dec, + pointing_icrs.icrs.ra, + pointing_icrs.icrs.dec, ) else: raise ValueError( From 4519509cb77458f3e6de638e6ceced619abc80a9 Mon Sep 17 00:00:00 2001 From: Quentin Remy Date: Wed, 17 Jul 2024 17:59:06 +0200 Subject: [PATCH 06/55] change sign to match previous behaviour ??? Signed-off-by: Quentin Remy --- gammapy/makers/utils.py | 1 + 1 file changed, 1 insertion(+) diff --git a/gammapy/makers/utils.py b/gammapy/makers/utils.py index 17b3476a43..f78964ced0 100644 --- a/gammapy/makers/utils.py +++ b/gammapy/makers/utils.py @@ -89,6 +89,7 @@ def _get_fov_coords(pointing, irf, geom, use_region_center=True, obstime=None): pointing_icrs.icrs.ra, pointing_icrs.icrs.dec, ) + fov_lon = -fov_lon else: raise ValueError( f"Unsupported background coordinate system: {irf.fov_alignment!r}" From c860563767e00136df96ca35e77949d1d9e5d311 Mon Sep 17 00:00:00 2001 From: Quentin Remy Date: Fri, 19 Jul 2024 19:05:50 +0200 Subject: [PATCH 07/55] add REVERSE_FOV_LON global Signed-off-by: Quentin Remy --- gammapy/makers/utils.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/gammapy/makers/utils.py b/gammapy/makers/utils.py index f78964ced0..0a7fea745a 100644 --- a/gammapy/makers/utils.py +++ b/gammapy/makers/utils.py @@ -14,6 +14,8 @@ from gammapy.utils.coordinates import sky_to_fov from gammapy.utils.regions import compound_region_to_regions +REVERSE_FOV_LON = False + __all__ = [ "make_counts_rad_max", "make_edisp_kernel_map", @@ -89,7 +91,8 @@ def _get_fov_coords(pointing, irf, geom, use_region_center=True, obstime=None): pointing_icrs.icrs.ra, pointing_icrs.icrs.dec, ) - fov_lon = -fov_lon + if REVERSE_FOV_LON: + fov_lon = -fov_lon else: raise ValueError( f"Unsupported background coordinate system: {irf.fov_alignment!r}" From 672fa9b8928d1ce0858d6d975a2d44c0c46f58b1 Mon Sep 17 00:00:00 2001 From: Quentin Remy Date: Tue, 30 Jul 2024 16:39:33 +0200 Subject: [PATCH 08/55] change FoVAlignment based on meta infos Signed-off-by: Quentin Remy --- gammapy/data/metadata.py | 2 ++ gammapy/data/observations.py | 24 ++++++++++++++++++++++-- gammapy/irf/core.py | 3 +++ gammapy/makers/tests/test_spectrum.py | 6 ++++++ gammapy/makers/utils.py | 6 ++---- 5 files changed, 35 insertions(+), 6 deletions(-) diff --git a/gammapy/data/metadata.py b/gammapy/data/metadata.py index ae2d55f78e..c5033e4393 100644 --- a/gammapy/data/metadata.py +++ b/gammapy/data/metadata.py @@ -107,6 +107,8 @@ def from_header(cls, header, format="gadf"): "PRESSURE", "RELHUM", "NSBLEVEL", + "CREATOR", + "HDUVERS", ] optional = dict() for key in optional_keywords: diff --git a/gammapy/data/observations.py b/gammapy/data/observations.py index 953c08ba61..d51f18d71f 100644 --- a/gammapy/data/observations.py +++ b/gammapy/data/observations.py @@ -16,6 +16,7 @@ from astropy.utils import lazyproperty import matplotlib.pyplot as plt from gammapy.utils.deprecation import GammapyDeprecationWarning +from gammapy.irf import FoVAlignment from gammapy.utils.fits import LazyFitsData, earth_location_to_dict from gammapy.utils.metadata import CreatorMetaData, TargetMetaData, TimeInfoMetaData from gammapy.utils.scripts import make_path @@ -65,7 +66,7 @@ class Observation: aeff = LazyFitsData(cache=False) edisp = LazyFitsData(cache=False) psf = LazyFitsData(cache=False) - bkg = LazyFitsData(cache=False) + _bkg = LazyFitsData(cache=False) _rad_max = LazyFitsData(cache=True) _events = LazyFitsData(cache=False) _gti = LazyFitsData(cache=True) @@ -90,7 +91,7 @@ def __init__( self.aeff = aeff self.edisp = edisp self.psf = psf - self.bkg = bkg + self._bkg = bkg self._rad_max = rad_max self._gti = gti self._events = events @@ -105,6 +106,25 @@ def _repr_html_(self): except AttributeError: return f"
{html.escape(str(self))}
" + @property + def bkg(self): + """Background of the observation.""" + bkg = self._bkg + # used for backward compatibility of old HESS data + if ( + bkg + and self._meta + and self._meta.optional + and self._meta.optional["CREATOR"] == "SASH FITS::EventListWriter" + and self.meta.optional["HDUVERS"] == "0.2" + ): + bkg._fov_alignment = FoVAlignment.REVERSE_LON_RADEC + return bkg + + @bkg.setter + def bkg(self, value): + self._bkg = value + @property def meta(self): """Return metadata container.""" diff --git a/gammapy/irf/core.py b/gammapy/irf/core.py index 3aa483af6b..d0f031ad71 100644 --- a/gammapy/irf/core.py +++ b/gammapy/irf/core.py @@ -33,6 +33,9 @@ class FoVAlignment(str, Enum): ALTAZ = "ALTAZ" RADEC = "RADEC" + REVERSE_LON_RADEC = ( + "REVERSE_LON_RADEC" # used for backward compatibility of old HESS data + ) class IRF(metaclass=abc.ABCMeta): diff --git a/gammapy/makers/tests/test_spectrum.py b/gammapy/makers/tests/test_spectrum.py index ae46443331..4949fcb24a 100644 --- a/gammapy/makers/tests/test_spectrum.py +++ b/gammapy/makers/tests/test_spectrum.py @@ -7,6 +7,7 @@ from regions import CircleSkyRegion from gammapy.data import DataStore from gammapy.datasets import SpectrumDataset +from gammapy.irf import FoVAlignment from gammapy.makers import ( ReflectedRegionsBackgroundMaker, SafeMaskMaker, @@ -125,6 +126,11 @@ def test_spectrum_dataset_maker_hess_dl3(spectrum_dataset_crab, observations_hes datasets = [] for obs in observations_hess_dl3: + + assert obs.meta.optional["CREATOR"] == "SASH FITS::EventListWriter" + assert obs.meta.optional["HDUVERS"] == "0.2" + assert obs.bkg.fov_alignment == FoVAlignment.REVERSE_LON_RADEC + dataset = maker.run(spectrum_dataset_crab, obs) datasets.append(dataset) diff --git a/gammapy/makers/utils.py b/gammapy/makers/utils.py index 0a7fea745a..da25161118 100644 --- a/gammapy/makers/utils.py +++ b/gammapy/makers/utils.py @@ -14,8 +14,6 @@ from gammapy.utils.coordinates import sky_to_fov from gammapy.utils.regions import compound_region_to_regions -REVERSE_FOV_LON = False - __all__ = [ "make_counts_rad_max", "make_edisp_kernel_map", @@ -84,14 +82,14 @@ def _get_fov_coords(pointing, irf, geom, use_region_center=True, obstime=None): fov_lon, fov_lat = sky_to_fov( altaz_coord.az, altaz_coord.alt, pointing_altaz.az, pointing_altaz.alt ) - elif irf.fov_alignment == FoVAlignment.RADEC: + elif irf.fov_alignment in [FoVAlignment.RADEC, FoVAlignment.REVERSE_LON_RADEC]: fov_lon, fov_lat = sky_to_fov( sky_coord.icrs.ra, sky_coord.icrs.dec, pointing_icrs.icrs.ra, pointing_icrs.icrs.dec, ) - if REVERSE_FOV_LON: + if irf.fov_alignment == FoVAlignment.REVERSE_LON_RADEC: fov_lon = -fov_lon else: raise ValueError( From cf9a99917c8169643f7fcd64898e07ac90e52213 Mon Sep 17 00:00:00 2001 From: Quentin Remy Date: Tue, 30 Jul 2024 16:57:06 +0200 Subject: [PATCH 09/55] fix Signed-off-by: Quentin Remy --- gammapy/data/observations.py | 21 ++++++++++++--------- 1 file changed, 12 insertions(+), 9 deletions(-) diff --git a/gammapy/data/observations.py b/gammapy/data/observations.py index d51f18d71f..fbd2b8fd80 100644 --- a/gammapy/data/observations.py +++ b/gammapy/data/observations.py @@ -111,14 +111,17 @@ def bkg(self): """Background of the observation.""" bkg = self._bkg # used for backward compatibility of old HESS data - if ( - bkg - and self._meta - and self._meta.optional - and self._meta.optional["CREATOR"] == "SASH FITS::EventListWriter" - and self.meta.optional["HDUVERS"] == "0.2" - ): - bkg._fov_alignment = FoVAlignment.REVERSE_LON_RADEC + try: + if ( + bkg + and self._meta + and self._meta.optional + and self._meta.optional["CREATOR"] == "SASH FITS::EventListWriter" + and self._meta.optional["HDUVERS"] == "0.2" + ): + bkg._fov_alignment = FoVAlignment.REVERSE_LON_RADEC + except KeyError: + pass return bkg @bkg.setter @@ -157,7 +160,7 @@ def rad_max(self): def available_hdus(self): """Which HDUs are available.""" available_hdus = [] - keys = ["_events", "_gti", "aeff", "edisp", "psf", "bkg", "_rad_max"] + keys = ["_events", "_gti", "aeff", "edisp", "psf", "_bkg", "_rad_max"] hdus = ["events", "gti", "aeff", "edisp", "psf", "bkg", "rad_max"] for key, hdu in zip(keys, hdus): available = self.__dict__.get(key, False) From 224de8008ed8a81df53cca9245e8a95aff3f30be Mon Sep 17 00:00:00 2001 From: Quentin Remy Date: Tue, 30 Jul 2024 17:06:20 +0200 Subject: [PATCH 10/55] Update gammapy/irf/core.py Co-authored-by: Maximilian Linhoff --- gammapy/irf/core.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/gammapy/irf/core.py b/gammapy/irf/core.py index d0f031ad71..91d7ddf52d 100644 --- a/gammapy/irf/core.py +++ b/gammapy/irf/core.py @@ -33,9 +33,8 @@ class FoVAlignment(str, Enum): ALTAZ = "ALTAZ" RADEC = "RADEC" - REVERSE_LON_RADEC = ( - "REVERSE_LON_RADEC" # used for backward compatibility of old HESS data - ) + # used for backward compatibility of old HESS data + REVERSE_LON_RADEC = "REVERSE_LON_RADEC" class IRF(metaclass=abc.ABCMeta): From f3fde440b5a5c7111d78b48ca41ab8d62461226c Mon Sep 17 00:00:00 2001 From: Astro-Kirsty Date: Mon, 21 Oct 2024 16:11:04 +0200 Subject: [PATCH 11/55] Link to the safe maker api Signed-off-by: Astro-Kirsty --- examples/tutorials/api/makers.py | 25 +++++++++++++------------ gammapy/makers/safe.py | 6 ++++-- 2 files changed, 17 insertions(+), 14 deletions(-) diff --git a/examples/tutorials/api/makers.py b/examples/tutorials/api/makers.py index 866517b99e..e13bfd0564 100644 --- a/examples/tutorials/api/makers.py +++ b/examples/tutorials/api/makers.py @@ -117,12 +117,13 @@ # `background_oversampling` parameter that defines the oversampling # factor in energy used to compute the background (default is None). # +# .. _safe-data-range: # Safe data range handling # ------------------------ # -# To exclude the data range from a `~gammapy.makers.MapDataset`, that is associated with -# high systematics on instrument response functions, a `~gammapy.makers.MapDataset.mask_safe` -# can be defined. The `~gammapy.makers.MapDataset.mask_safe` is a `~gammapy.maps.Map` object +# To exclude the data range from a `~gammapy.datasets.MapDataset`, that is associated with +# high systematics on instrument response functions, a `~gammapy.datasets.MapDataset.mask_safe` +# can be defined. The `~gammapy.datasets.MapDataset.mask_safe` is a `~gammapy.maps.Map` object # with `bool` data type, which indicates for each pixel, whether it should be included in # the analysis. The convention is that a value of `True` or `1` # includes the pixel, while a value of `False` or `0` excludes a @@ -141,8 +142,8 @@ # observation center are excluded # - bkg-peak, the energy threshold is defined as the upper edge of the # energy bin with the highest predicted background rate. This method -# was introduced in the HESS DL3 validation paper: -# https://arxiv.org/pdf/1910.08088.pdf +# was introduced in the +# `H.E.S.S. DL3 validation paper `__ # # Note that currently some methods computing a safe energy range # ("aeff-default", "aeff-max" and "edisp-bias") determine a true energy range and @@ -175,7 +176,7 @@ # Background estimation # --------------------- # -# The background computed by the `MapDatasetMaker` gives the number of +# The background computed by the `~gammapy.makers.MapDatasetMaker` gives the number of # counts predicted by the background IRF of the observation. Because its # actual normalization, or even its spectral shape, might be poorly # constrained, it is necessary to correct it with the data themselves. @@ -241,7 +242,7 @@ # ------------------- # # The data reduction steps can be combined in a single loop to run a full -# data reduction chain. For this the `MapDatasetMaker` is run first and +# data reduction chain. For this the `~gammapy.makers.MapDatasetMaker` is run first and # the output dataset is the passed on to the next maker step. Finally, the # dataset per observation is stacked into a larger map. # @@ -273,7 +274,7 @@ ###################################################################### # To maintain good performance it is always recommended to do a cutout of -# the `MapDataset` as shown above. In case you want to increase the +# the `~gammapy.datasets.MapDataset` as shown above. In case you want to increase the # offset-cut later, you can also choose a larger width of the cutout than # `2 * offset_max`. # @@ -308,8 +309,8 @@ # count spectra, background is implicitly modeled via the OFF counts # spectrum. # -# The `~gammapy.datasets.SpectrumDatasetMaker` make spectrum dataset for a single -# observation. In that case the irfs and background are computed at a +# The `~gammapy.makers.SpectrumDatasetMaker` make spectrum dataset for a single +# observation. In that case the IRFs and background are computed at a # single fixed offset, which is recommended only for point-sources. # # Here is an example of data reduction loop to create @@ -349,7 +350,7 @@ # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ # It can often be useful to know the total number of hours spent # in the given field of view (without correcting for the acceptance -# variation). This can be computed using `make_observation_time_map` +# variation). This can be computed using `~gammapy.makers.utils.make_observation_time_map` # as shown below # @@ -420,7 +421,7 @@ ###################################################################### # To get the value of the observation time at a particular position, -# use `get_by_coord` +# use ``get_by_coord`` obs_time_src = total_obstime.get_by_coord(source_pos) effective_times_src = effective_livetime.get_by_coord( diff --git a/gammapy/makers/safe.py b/gammapy/makers/safe.py index 73d9b60b8c..54e1635500 100644 --- a/gammapy/makers/safe.py +++ b/gammapy/makers/safe.py @@ -17,9 +17,11 @@ class SafeMaskMaker(Maker): """Make safe data range mask for a given observation. + For more information see :ref:`safe-data-range`. + .. warning:: - Currently some methods computing a safe energy range ("aeff-default", + Currently, some methods computing a safe energy range ("aeff-default", "aeff-max" and "edisp-bias") determine a true energy range and apply it to reconstructed energy, effectively neglecting the energy dispersion. @@ -280,7 +282,7 @@ def make_mask_energy_bkg_peak(self, dataset, observation=None): The energy threshold is defined as the lower edge of the energy bin with the highest predicted background rate. This is to ensure analysis in a region where a Powerlaw approximation to the background spectrum is valid. - The is motivated by its use in the HESS DL3 + The is motivated by its use in the H.E.S.S. DL3 validation paper: https://arxiv.org/pdf/1910.08088.pdf Parameters From 91c89f0a4630d745c81d06b0f4b3b7b7f9ee2f39 Mon Sep 17 00:00:00 2001 From: Daniel Morcuende Date: Mon, 21 Oct 2024 19:01:23 +0200 Subject: [PATCH 12/55] remove light/dark theme switch from docs sphinx config --- docs/conf.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/conf.py b/docs/conf.py index 86a2edbc91..cf3ee50360 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -237,7 +237,7 @@ def setup(app): "json_url": "https://docs.gammapy.org/stable/switcher.json", "version_match": switch_version, }, - "navbar_end": ["version-switcher", "theme-switcher", "navbar-icon-links"], + "navbar_end": ["version-switcher", "navbar-icon-links"], "navigation_with_keys": True, # footers "footer_start": ["copyright"], From 6bf10130eb89aff707895cc2519b09f0fec02197 Mon Sep 17 00:00:00 2001 From: Astro-Kirsty Date: Tue, 15 Oct 2024 14:46:34 +0200 Subject: [PATCH 13/55] Fix docs build of the examples in estimators rst Signed-off-by: Astro-Kirsty --- docs/user-guide/estimators.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/user-guide/estimators.rst b/docs/user-guide/estimators.rst index 76bc6843d0..4cbe4db47e 100644 --- a/docs/user-guide/estimators.rst +++ b/docs/user-guide/estimators.rst @@ -91,7 +91,7 @@ arguments are translated into a TS difference assuming ``ts = n_sigma ** 2``. .. _sedtypes: SED types -~~~~~~~~~ +^^^^^^^^^ In addition to the norm values a reference spectral model and energy ranges are given. Using this reference spectral model the norm values can be converted From 99828903424d7b18eff0d09828f643e345d737bc Mon Sep 17 00:00:00 2001 From: Astro-Kirsty Date: Wed, 16 Oct 2024 12:10:52 +0200 Subject: [PATCH 14/55] Adjust tute Signed-off-by: Astro-Kirsty --- examples/tutorials/data/hess.py | 1 + 1 file changed, 1 insertion(+) diff --git a/examples/tutorials/data/hess.py b/examples/tutorials/data/hess.py index 8337b12e3f..0c6131715b 100644 --- a/examples/tutorials/data/hess.py +++ b/examples/tutorials/data/hess.py @@ -1,6 +1,7 @@ """ H.E.S.S. with Gammapy ===================== + Explore H.E.S.S. event lists and IRFs. From c64f28102a46fd5e4e351e94b397a05717a43d1c Mon Sep 17 00:00:00 2001 From: Astro-Kirsty Date: Fri, 18 Oct 2024 15:32:36 +0200 Subject: [PATCH 15/55] fix links in rsts Signed-off-by: Astro-Kirsty --- docs/user-guide/astro/darkmatter/index.rst | 8 ++--- docs/user-guide/astro/index.rst | 4 +-- docs/user-guide/catalog.rst | 26 +++++++-------- docs/user-guide/datasets/index.rst | 18 +++++----- docs/user-guide/dl3.rst | 4 +-- docs/user-guide/estimators.rst | 6 ++-- docs/user-guide/hli.rst | 2 +- docs/user-guide/makers/index.rst | 2 +- docs/user-guide/maps/index.rst | 38 +++++++++++----------- docs/user-guide/scripts/index.rst | 2 +- docs/user-guide/utils.rst | 6 ++-- docs/user-guide/visualization/index.rst | 2 +- 12 files changed, 60 insertions(+), 58 deletions(-) diff --git a/docs/user-guide/astro/darkmatter/index.rst b/docs/user-guide/astro/darkmatter/index.rst index daf81c33cb..13808b86d9 100644 --- a/docs/user-guide/astro/darkmatter/index.rst +++ b/docs/user-guide/astro/darkmatter/index.rst @@ -89,11 +89,11 @@ module `DarkBit`_. DarkBit can be used for directly computing observables and likelihoods, for any combination of parameter values in some underlying particle model. -Gammapy tutorial -================ +Using gammapy.astro.darkmatter +------------------------------ + +.. minigallery:: gammapy.astro.darkmatter -.. minigallery:: `gammapy.astro.darkmatter.DarkMatterAnnihilationSpectralModel` - :add-heading: .. _Cirelli et al. 2011: http://iopscience.iop.org/article/10.1088/1475-7516/2011/03/051/pdf diff --git a/docs/user-guide/astro/index.rst b/docs/user-guide/astro/index.rst index 9fcd287eed..3bc8c7184e 100644 --- a/docs/user-guide/astro/index.rst +++ b/docs/user-guide/astro/index.rst @@ -9,14 +9,14 @@ This module contains utility functions for some astrophysical scenarios: * `~gammapy.astro.population` for astrophysical population models * `~gammapy.astro.darkmatter` for dark matter spatial and spectral models -The ``gammapy.astro`` sub-package is in a prototyping phase and its scope and future +The `gammapy.astro` sub-package is in a prototyping phase and its scope and future are currently being discussed. It is likely that some functionality will be removed or split out into a separate package at some point. Getting started --------------- -The ``gammapy.astro`` namespace is empty. Use these import statements: +The `gammapy.astro` namespace is empty. Use these import statements: .. testcode:: diff --git a/docs/user-guide/catalog.rst b/docs/user-guide/catalog.rst index 6a4225629e..0a426671d7 100644 --- a/docs/user-guide/catalog.rst +++ b/docs/user-guide/catalog.rst @@ -5,17 +5,17 @@ Source catalogs `gammapy.catalog` provides convenient access to common gamma-ray source catalogs. -* ``hgps`` / `SourceCatalogHGPS` - H.E.S.S. Galactic plane survey (HGPS) -* ``gamma-cat`` / `SourceCatalogGammaCat` - An open catalog of gamma-ray sources -* ``3fgl`` / `SourceCatalog3FGL` - LAT 4-year point source catalog -* ``4fgl`` / `SourceCatalog4FGL` - LAT 8-year point source catalog -* ``2fhl`` / `SourceCatalog2FHL` - LAT second high-energy source catalog -* ``3fhl`` / `SourceCatalog3FHL` - LAT third high-energy source catalog -* ``2hwc`` / `SourceCatalog2HWC` - 2HWC catalog from the HAWC observatory -* ``3hwc`` / `SourceCatalog3HWC` - 3HWC catalog from the HAWC observatory - -For each catalog, a `SourceCatalog` class is provided to represent the catalog table, -and a matching `SourceCatalogObject` class to represent one catalog source and table row. +* ``hgps`` / `~gammapy.catalog.SourceCatalogHGPS` - H.E.S.S. Galactic plane survey (HGPS) +* ``gamma-cat`` / `~gammapy.catalog.SourceCatalogGammaCat` - An open catalog of gamma-ray sources +* ``3fgl`` / `~gammapy.catalog.SourceCatalog3FGL` - LAT 4-year point source catalog +* ``4fgl`` / `~gammapy.catalog.SourceCatalog4FGL` - LAT 8-year point source catalog +* ``2fhl`` / `~gammapy.catalog.SourceCatalog2FHL` - LAT second high-energy source catalog +* ``3fhl`` / `~gammapy.catalog.SourceCatalog3FHL` - LAT third high-energy source catalog +* ``2hwc`` / `~gammapy.catalog.SourceCatalog2HWC` - 2HWC catalog from the HAWC observatory +* ``3hwc`` / `~gammapy.catalog.SourceCatalog3HWC` - 3HWC catalog from the HAWC observatory + +For each catalog, a `~gammapy.catalog.SourceCatalog` class is provided to represent the catalog table, +and a matching `~gammapy.catalog.SourceCatalogObject` class to represent one catalog source and table row. The main functionality provided is methods that map catalog information to `~gammapy.modeling.models.SkyModel`, `~gammapy.modeling.models.SpectralModel`, @@ -29,9 +29,9 @@ or to create initial source models for certain energy bands and sky regions. Using gammapy.catalog --------------------- -.. minigallery:: gammapy.catalog.SourceCatalog3FHL +.. minigallery:: `~gammapy.catalog.SourceCatalog3FHL` :add-heading: -.. minigallery:: gammapy.catalog.SourceCatalogGammaCat +.. minigallery:: `~gammapy.catalog.SourceCatalogGammaCat` :add-heading: diff --git a/docs/user-guide/datasets/index.rst b/docs/user-guide/datasets/index.rst index 19f73627b0..2bd4fdf715 100644 --- a/docs/user-guide/datasets/index.rst +++ b/docs/user-guide/datasets/index.rst @@ -8,7 +8,7 @@ Datasets (DL4) The `gammapy.datasets` sub-package contains classes to handle reduced gamma-ray data for modeling and fitting. -The `Dataset` class bundles reduced data, IRFs and model to perform +The `~gammapy.datasets.Dataset` class bundles reduced data, IRFs and model to perform likelihood fitting and joint-likelihood fitting. All datasets contain a `~gammapy.modeling.models.Models` container with one or more `~gammapy.modeling.models.SkyModel` objects that represent additive emission @@ -96,7 +96,7 @@ a joint fit across multiple datasets. Predicted counts ---------------- -The total number of predicted counts from a `MapDataset` are computed per bin like: +The total number of predicted counts from a `~gammapy.datasets.MapDataset` are computed per bin like: .. math:: @@ -174,7 +174,7 @@ Here, :math:`k` denotes a bin in reconstructed energy, For the model evaluation, an important factor that needs to be accounted for is that the energy threshold changes between observations. -With the above implementation using a `~gammapy.irf.EDispersionMap`, +With the above implementation using a `~gammapy.irf.EDispMap`, the `npred` is conserved, ie, the predicted number of counts on the stacked dataset is the sum expected by stacking the `npred` of the individual runs, @@ -220,11 +220,13 @@ stack runs taken under similar conditions and then do a joint fit on the stacked Serialisation of datasets ------------------------- -The various `Map` objects contained in `~gammapy.datasets.MapDataset` and `~gammapy.datasets.MapDatasetOnOff` are serialised according to `GADF Sky Maps `__. +The various `~gammapy.maps.Map` objects contained in `~gammapy.datasets.MapDataset` and +`~gammapy.datasets.MapDatasetOnOff` are serialised according to +`GADF Sky Maps `__. A hdulist is created with the different attributes, and each of these are written with the data -contained in a `BinTableHDU` with a `WcsGeom` and a `BANDS HDU` specifying the non-spatial dimensions. -Optionally, a `meta_table` is also written as an `astropy.table.Table` containing various information -about the observations which created the dataset. While the `meta_table` can contain useful information for +contained in a ``BinTableHDU`` with a `~gammapy.maps.WcsGeom` and a ``BANDS HDU`` specifying the non-spatial dimensions. +Optionally, a ``meta_table`` is also written as an `astropy.table.Table` containing various information +about the observations which created the dataset. While the ``meta_table`` can contain useful information for later stage analysis, it is not used anywhere internally within gammapy. `~gammapy.datasets.SpectrumDataset` follows a similar convention as for `~gammapy.datasets.MapDataset`, but uses a @@ -233,7 +235,7 @@ later stage analysis, it is not used anywhere internally within gammapy. either according to the above specification, or (by default), according to the `OGIP standards `__. -`~gammapy.datasets.FluxPointsDatasets` are serialised as `gammapy.estimators.FluxPoints` objects, which contains +`~gammapy.datasets.FluxPointsDataset` are serialised as `gammapy.estimators.FluxPoints` objects, which contains a set of `gammapy.maps.Map` objects storing the estimated flux as function of energy, and some optional quantities like typically errors, upper limits, etc. It also contains a reference model, serialised as a `~gammapy.modeling.models.TemplateSpectralModel`. diff --git a/docs/user-guide/dl3.rst b/docs/user-guide/dl3.rst index abea89a040..9eb8cb9dbc 100644 --- a/docs/user-guide/dl3.rst +++ b/docs/user-guide/dl3.rst @@ -32,7 +32,7 @@ You can use the `~gammapy.data.EventList` class to load IACT gamma-ray event lis filename = '$GAMMAPY_DATA/hess-dl3-dr1/data/hess_dl3_dr1_obs_id_023523.fits.gz' events = EventList.read(filename) -To load Fermi-LAT event lists, use the `~gammapy.data.EventListLAT` class: +To load Fermi-LAT event lists, you can also use the `~gammapy.data.EventList` class: .. testcode:: @@ -139,7 +139,7 @@ Writing event lists and GTIs to file To write the events or GTIs separately, one can just save the underlying `astropy.table.Table`. To have an event file written in a correct DL3 format, it is -necessary to utilise the ``write`` method available for`~gammapy.data.Observation`. +necessary to utilise the ``write`` method available for `~gammapy.data.Observation`. It will write the `~gammapy.data.EventList` and their associated GTIs together in the same FITS file according to the format specifications. To avoid writing IRFs along the ``EventList`` one has to set ``include_irfs`` to ``False``. See the example below: diff --git a/docs/user-guide/estimators.rst b/docs/user-guide/estimators.rst index 4cbe4db47e..5c87e5f901 100644 --- a/docs/user-guide/estimators.rst +++ b/docs/user-guide/estimators.rst @@ -35,7 +35,7 @@ In general the flux can be estimated using two methods: not re-optimising other parameters, one can estimate the significance based on the analytical solution by [LiMa1983]. In this case the "best fit" flux and significance are given by the excess over the null hypothesis. This method is also named - **backward folding**. This method is currently only exposed in the `ExcessMapEstimator` + **backward folding**. This method is currently only exposed in the `~gammapy.estimators.ExcessMapEstimator` Energy edges @@ -112,7 +112,7 @@ e2dnde Differential energy flux at ``e_ref`` The same can be applied for the error and upper limit information. More information can be found on the `likelihood SED type page`_. -The `FluxPoints` and `FluxMaps` objects can optionally define meta +The `~gammapy.estimators.FluxPoints` and `~gammapy.estimators.FluxMaps` objects can optionally define meta data with the following valid keywords: ================= ================================================= @@ -141,7 +141,7 @@ A note on negative flux and upper limit values: Flux maps --------- -This how to compute flux maps with the `ExcessMapEstimator`: +This how to compute flux maps with the `~gammapy.estimators.ExcessMapEstimator`: .. testcode:: diff --git a/docs/user-guide/hli.rst b/docs/user-guide/hli.rst index ad0c8a7e06..1fcd585da5 100644 --- a/docs/user-guide/hli.rst +++ b/docs/user-guide/hli.rst @@ -14,7 +14,7 @@ the functionality that is available in the sub-packages to support standard analysis use case in a convenient way. -Using gammapy.datasets +Using gammapy.analysis ---------------------- .. minigallery:: gammapy.analysis.Analysis diff --git a/docs/user-guide/makers/index.rst b/docs/user-guide/makers/index.rst index e67091add6..1a8fe24b72 100644 --- a/docs/user-guide/makers/index.rst +++ b/docs/user-guide/makers/index.rst @@ -24,7 +24,7 @@ Background estimation Safe data range definition -------------------------- -The definition of a safe data range is done using the `SafeMaskMaker` or manually. +The definition of a safe data range is done using the `~gammapy.makers.SafeMaskMaker` or manually. Using gammapy.makers diff --git a/docs/user-guide/maps/index.rst b/docs/user-guide/maps/index.rst index 9b783254e8..d272251eea 100644 --- a/docs/user-guide/maps/index.rst +++ b/docs/user-guide/maps/index.rst @@ -20,11 +20,11 @@ non-spatial dimensions and can represent images (2D), cubes (3D), or hypercubes `gammapy.maps` is organized around two data structures: *geometry* classes -inheriting from `~Geom` and *map* classes inheriting from `~Map`. A geometry -defines the map boundaries, pixelization scheme, and provides methods for -converting to/from map and pixel coordinates. A map owns a `~Geom` instance +inheriting from `~gammapy.maps.Geom` and *map* classes inheriting from `~gammapy.maps.Map`. A geometry +defines the map boundaries, pixelisation scheme, and provides methods for +converting to/from map and pixel coordinates. A map owns a `~gammapy.maps.Geom` instance as well as a data array containing map values. Where possible it is recommended -to use the abstract `~Map` interface for accessing or updating the contents of a +to use the abstract `~gammapy.maps.Map` interface for accessing or updating the contents of a map as this allows algorithms to be used interchangeably with different map representations. The following reviews methods of the abstract map interface. @@ -33,11 +33,11 @@ Getting started with maps ------------------------- All map objects have an abstract interface provided through the methods of the -`~Map`. These methods can be used for accessing and manipulating the contents of +`~gammapy.maps.Map`. These methods can be used for accessing and manipulating the contents of a map without reference to the underlying data representation (e.g. whether a -map uses WCS or HEALPix pixelization). For applications which do depend on the +map uses WCS or HEALPix pixelisation). For applications which do depend on the specific representation one can also work directly with the classes derived from -`~Map`. In the following we review some of the basic methods for working with +`~gammapy.maps.Map`. In the following we review some of the basic methods for working with map objects, more details are given in the :doc:`/tutorials/api/maps` tutorial. @@ -47,20 +47,20 @@ Accessor methods ---------------- All map objects have a set of accessor methods provided through the abstract -`~Map` class. These methods can be used to access or update the contents of the +`~gammapy.maps.Map` class. These methods can be used to access or update the contents of the map irrespective of its underlying representation. Four types of accessor methods are provided: * ``get`` : Return the value of the map at the pixel containing the - given coordinate (`~Map.get_by_idx`, `~Map.get_by_pix`, `~Map.get_by_coord`). + given coordinate (`~gammapy.maps.Map.get_by_idx`, `~gammapy.maps.Map.get_by_pix`, `~gammapy.maps.Map.get_by_coord`). * ``interp`` : Interpolate or extrapolate the value of the map at an arbitrary coordinate. * ``set`` : Set the value of the map at the pixel containing the - given coordinate (`~Map.set_by_idx`, `~Map.set_by_pix`, `~Map.set_by_coord`). + given coordinate (`~gammapy.maps.Map.set_by_idx`, `~gammapy.maps.Map.set_by_pix`, `~gammapy.maps.Map.set_by_coord`). * ``fill`` : Increment the value of the map at the pixel containing the given coordinate with a unit weight or the value in the optional - ``weights`` argument (`~Map.fill_by_idx`, `~Map.fill_by_pix`, - `~Map.fill_by_coord`). + ``weights`` argument (`~gammapy.maps.Map.fill_by_idx`, `~gammapy.maps.Map.fill_by_pix`, + `~gammapy.maps.Map.fill_by_coord`). Accessor methods accept as their first argument a coordinate tuple containing scalars, lists, or numpy arrays with one tuple element for each dimension of the @@ -84,7 +84,7 @@ coordinates can be expressed in one of three coordinate systems: array for each non-spatial dimension. The coordinate system accepted by a given accessor method can be inferred from -the suffix of the method name (e.g. `~Map.get_by_idx`). The following +the suffix of the method name (e.g. `~gammapy.maps.Map.get_by_idx`). The following demonstrates how one can access the same pixels of a WCS map using each of the three coordinate systems: @@ -141,7 +141,7 @@ following demonstrates how one can set pixel values: Interface with MapCoord and SkyCoord ------------------------------------ -The ``coord`` accessor methods accept `dict`, `~gammapy.maps.MapCoord`, and +The ``coord`` accessor methods accept ``dict``, `~gammapy.maps.MapCoord`, and `~astropy.coordinates.SkyCoord` arguments in addition to the standard `tuple` of `~numpy.ndarray` argument. When using a `tuple` argument a `~astropy.coordinates.SkyCoord` can be used instead of longitude and latitude @@ -166,7 +166,7 @@ transformed to match the coordinate system of the map. m.set_by_coord((skycoord, energy), [0.5, 1.5]) m.get_by_coord((skycoord, energy)) -A `~gammapy.maps.MapCoord` or `dict` argument can be used to interact with a map object +A `~gammapy.maps.MapCoord` or ``dict`` argument can be used to interact with a map object without reference to the axis ordering of the map geometry: .. testcode:: @@ -188,9 +188,9 @@ MapCoord `~gammapy.maps.MapCoord` is an N-dimensional coordinate object that stores both spatial and non-spatial coordinates and is accepted by all ``coord`` methods. A `~gammapy.maps.MapCoord` -can be created with or without explicitly named axes with `MapCoord.create`. +can be created with or without explicitly named axes with `~gammapy.maps.MapCoord.create`. Axes of a `~gammapy.maps.MapCoord` can be accessed by index, name, or attribute. A `~gammapy.maps.MapCoord` -without explicit axis names can be created by calling `MapCoord.create` with a +without explicit axis names can be created by calling `~gammapy.maps.MapCoord.create` with a `tuple` argument: .. testcode:: @@ -231,8 +231,8 @@ latitude. Non-spatial axes are assigned a default name ``axis{I}`` where ``{I}`` is the index of the non-spatial dimension. `~gammapy.maps.MapCoord` objects created without named axes must have the same axis ordering as the map geometry. -A `~gammapy.maps.MapCoord` with named axes can be created by calling `MapCoord.create` -with a `dict`: +A `~gammapy.maps.MapCoord` with named axes can be created by calling `~gammapy.maps.MapCoord.create` +with a ``dict``: .. testcode:: diff --git a/docs/user-guide/scripts/index.rst b/docs/user-guide/scripts/index.rst index 3ce4c0421e..63d299370e 100644 --- a/docs/user-guide/scripts/index.rst +++ b/docs/user-guide/scripts/index.rst @@ -150,7 +150,7 @@ Write your own CLI This section explains how to write your own command line interface (CLI). We will focus on the command line aspect, and use a very simple example where we -just call `gammapy.stats.CashCountsStatistics.sqrt_ts`. +just call `gammapy.stats.CashCountsStatistic.sqrt_ts`. From the interactive Python or IPython prompt or from a Jupyter notebook you just import the functionality you need and call it, like this: diff --git a/docs/user-guide/utils.rst b/docs/user-guide/utils.rst index dcb01e7dcd..b692ffcf7d 100644 --- a/docs/user-guide/utils.rst +++ b/docs/user-guide/utils.rst @@ -5,12 +5,12 @@ Utility functions ================= -``gammapy.utils`` is a collection of utility functions that are used in many +`gammapy.utils` is a collection of utility functions that are used in many places or don't fit in one of the other packages. -Since the various sub-modules of ``gammapy.utils`` are mostly unrelated, they +Since the various sub-modules of `gammapy.utils` are mostly unrelated, they are not imported into the top-level namespace. Here are some examples of how to -import functionality from the ``gammapy.utils`` sub-modules: +import functionality from the `gammapy.utils` sub-modules: .. testcode:: diff --git a/docs/user-guide/visualization/index.rst b/docs/user-guide/visualization/index.rst index 6b0f3a3b90..7583c205c3 100644 --- a/docs/user-guide/visualization/index.rst +++ b/docs/user-guide/visualization/index.rst @@ -19,7 +19,7 @@ used in gamma-ray astronomy (`colormap_hess` and `colormap_milagro`). Survey panel plots ------------------ -The `MapPanelPlotter` allows to split e.g. a galactic plane survey image with +The `~gammapy.visualization.MapPanelPlotter` allows to split e.g. a galactic plane survey image with a large aspect ratio into multiple panels. Here is a short example: .. plot:: From f6ceb1621954f8c62f54705fe2f59eb6a2864ba4 Mon Sep 17 00:00:00 2001 From: Maxime Regeard Date: Tue, 22 Oct 2024 12:38:29 +0200 Subject: [PATCH 16/55] correct docstring --- gammapy/datasets/flux_points.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/gammapy/datasets/flux_points.py b/gammapy/datasets/flux_points.py index 30e767527a..81915d0b91 100644 --- a/gammapy/datasets/flux_points.py +++ b/gammapy/datasets/flux_points.py @@ -40,7 +40,7 @@ def _get_reference_model(model, energy_bounds, margin_percent=70): class FluxPointsDataset(Dataset): - """Bundle a set of flux points with a parametric model, to compute fit statistic function using chi2 statistics. + """Bundle a set of flux points with a parametric model, to compute fit statistic function using different statistics (see ``stat_type``). For more information see :ref:`datasets`. From 8eb3067ec29b79e8d6147ec2bdedad76080c7c9d Mon Sep 17 00:00:00 2001 From: Atreyee Date: Tue, 22 Oct 2024 12:58:55 +0200 Subject: [PATCH 17/55] improve notebook Signed-off-by: Atreyee Signed-off-by: --- .../analysis-time/light_curve_flare.py | 57 +++++++++++-------- 1 file changed, 34 insertions(+), 23 deletions(-) diff --git a/examples/tutorials/analysis-time/light_curve_flare.py b/examples/tutorials/analysis-time/light_curve_flare.py index 573d9f17f3..679750c25e 100644 --- a/examples/tutorials/analysis-time/light_curve_flare.py +++ b/examples/tutorials/analysis-time/light_curve_flare.py @@ -86,6 +86,7 @@ ) from gammapy.maps import MapAxis, RegionGeom from gammapy.modeling.models import PowerLawSpectralModel, SkyModel +from gammapy.modeling import Fit ###################################################################### # Check setup @@ -126,9 +127,8 @@ # Define time intervals # --------------------- # -# We create the list of time intervals. Each time interval is an +# We create the list of time intervals, each of duration 10 mins. Each time interval is an # `astropy.time.Time` object, containing a start and stop time. -# t0 = Time("2006-07-29T20:30") duration = 10 * u.min @@ -236,26 +236,20 @@ ###################################################################### -# Define the Model +# Define underlying model # ---------------- # -# The actual flux will depend on the spectral shape assumed. For -# simplicity, we use the power law spectral model of index 3.4 used in the +# Since we use forward folding to obtain the flux points in each bin, exact values will depend on the underlying model. In this example, we use a power law as used in the # `reference -# paper `__. +# paper `. # -# Here we use only a spectral model in the -# `~gammapy.modeling.models.SkyModel` object. +# As we have are only using spectral datasats, we do not need any spatial models. # +# **Note** : All time bins must have the same spectral model. To see how to investigate spectral variability, see :doc:`time resolved spectroscopy notebook `. -spectral_model = PowerLawSpectralModel( - index=3.4, amplitude=2e-11 * u.Unit("1 / (cm2 s TeV)"), reference=1 * u.TeV -) -spectral_model.parameters["index"].frozen = False - +spectral_model = PowerLawSpectralModel(amplitude=1e-10 * u.Unit("1 / (cm2 s TeV)")) sky_model = SkyModel(spatial_model=None, spectral_model=spectral_model, name="pks2155") - ###################################################################### # Assign to model to all datasets # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -266,22 +260,29 @@ datasets.models = sky_model +# %%time +fit = Fit() +result = fit.run(datasets) +print(result.models.to_parameters_table()) + ###################################################################### # Extract the light curve # ----------------------- # # We first create the `~gammapy.estimators.LightCurveEstimator` for the # list of datasets we just produced. We give the estimator the name of the -# source component to be fitted. +# source component to be fitted. We can directly compute the light curve in multiple energy bins by supplying a list of `energy_edges`. +# # By default the likelihood scan is computed from 0.2 to 5.0. # Here, we increase the max value to 10.0, because we are # dealing with a large flare. lc_maker_1d = LightCurveEstimator( - energy_edges=[0.7, 20] * u.TeV, + energy_edges=[0.7, 1, 20] * u.TeV, source="pks2155", time_intervals=time_intervals, selection_optional="all", + n_jobs=4, ) lc_maker_1d.norm.scan_max = 10 @@ -301,21 +302,31 @@ # Finally we plot the result for the 1D lightcurve: # plt.figure(figsize=(8, 6)) -lc_1d.plot(marker="o") +lc_1d.plot(marker="o", axis_name="time", sed_type="flux") plt.show() ###################################################################### -# Light curves once obtained can be rebinned. -# Here, we rebin 4 adjacent bins together, to get 30 min bins +# Light curves once obtained can be rebinned using the likelihood profiles. +# Here, we rebin 3 adjacent bins together, to get 30 min bins. +# +# We will first slice `lc_1d` to obtain the lightcurve in the first energy bin # -axis_new = get_rebinned_axis(lc_1d, method="fixed-bins", group_size=3, axis_name="time") +slices = {"energy": slice(0, 1)} +sliced_lc = lc_1d.slice_by_idx(slices) +print(sliced_lc) + +axis_new = get_rebinned_axis( + sliced_lc, method="fixed-bins", group_size=3, axis_name="time" +) print(axis_new) -lc_new = lc_1d.resample_axis(axis_new) +lc_new = sliced_lc.resample_axis(axis_new) plt.figure(figsize=(8, 6)) -ax = lc_1d.plot(label="original") +ax = sliced_lc.plot(label="original") lc_new.plot(ax=ax, label="rebinned") plt.legend() -plt.show() \ No newline at end of file +plt.show() + +# We can use the sliced lightcurve to understand the variability, as shown in the doc:`/tutorials/analysis-time/variability_estimation` tutorial. From 6e5f790e48e6bbe1264e64841b874bada2143d22 Mon Sep 17 00:00:00 2001 From: Atreyee Date: Mon, 21 Oct 2024 16:18:30 +0200 Subject: [PATCH 18/55] improve ring notebook Signed-off-by: Atreyee Signed-off-by: --- .../tutorials/analysis-2d/ring_background.py | 58 ++++++++++--------- 1 file changed, 31 insertions(+), 27 deletions(-) diff --git a/examples/tutorials/analysis-2d/ring_background.py b/examples/tutorials/analysis-2d/ring_background.py index 7947fc923b..7ffb20e6df 100644 --- a/examples/tutorials/analysis-2d/ring_background.py +++ b/examples/tutorials/analysis-2d/ring_background.py @@ -39,7 +39,7 @@ """ -###################################################################### +# ##################################################################### # Setup # ----- # @@ -59,10 +59,11 @@ from gammapy.makers import RingBackgroundMaker from gammapy.visualization import plot_distribution from gammapy.utils.check import check_tutorials_setup + log = logging.getLogger(__name__) -###################################################################### +# ##################################################################### # Check setup # ----------- @@ -70,7 +71,7 @@ check_tutorials_setup() -###################################################################### +# ##################################################################### # Creating the config file # ------------------------ # @@ -104,11 +105,11 @@ config.datasets.geom.wcs.binsize = "0.02 deg" # Cutout size (for the run-wise event selection) -config.datasets.geom.selection.offset_max = 3.5 * u.deg +config.datasets.geom.selection.offset_max = 2.5 * u.deg # We now fix the energy axis for the counts map - (the reconstructed energy binning) config.datasets.geom.axes.energy.min = "0.5 TeV" -config.datasets.geom.axes.energy.max = "5 TeV" +config.datasets.geom.axes.energy.max = "10 TeV" config.datasets.geom.axes.energy.nbins = 10 # We need to extract the ring for each observation separately, hence, no stacking at this stage @@ -117,7 +118,7 @@ print(config) -###################################################################### +# ##################################################################### # Getting the reduced dataset # --------------------------- # @@ -144,7 +145,7 @@ analysis.get_datasets() -###################################################################### +# ##################################################################### # Extracting the ring background # ------------------------------ # @@ -154,7 +155,7 @@ # -###################################################################### +# ##################################################################### # Create exclusion mask # ~~~~~~~~~~~~~~~~~~~~~ # @@ -169,13 +170,13 @@ geom_image = geom.to_image().to_cube([energy_axis.squash()]) # Make the exclusion mask -regions = CircleSkyRegion(center=source_pos, radius=0.3 * u.deg) +regions = CircleSkyRegion(center=source_pos, radius=0.4 * u.deg) exclusion_mask = ~geom_image.region_mask([regions]) exclusion_mask.sum_over_axes().plot() plt.show() -###################################################################### +# ##################################################################### # For the present analysis, we use a ring with an inner radius of 0.5 deg # and width of 0.3 deg. # @@ -185,7 +186,7 @@ ) -###################################################################### +# ##################################################################### # Create a stacked dataset # ~~~~~~~~~~~~~~~~~~~~~~~~ # @@ -193,19 +194,19 @@ # together to create a single stacked map for further analysis # -# %%time energy_axis_true = analysis.datasets[0].exposure.geom.axes["energy_true"] stacked_on_off = MapDatasetOnOff.create( geom=geom_image, energy_axis_true=energy_axis_true, name="stacked" ) +# %%time for dataset in analysis.datasets: # Ring extracting makes sense only for 2D analysis dataset_on_off = ring_maker.run(dataset.to_image()) stacked_on_off.stack(dataset_on_off) -###################################################################### +# ##################################################################### # This `stacked_on_off` has `on` and `off` counts and acceptance # maps which we will use in all further analysis. The `acceptance` and # `acceptance_off` maps are the system acceptance of gamma-ray like @@ -215,7 +216,7 @@ print(stacked_on_off) -###################################################################### +# ##################################################################### # Compute correlated significance and correlated excess maps # ---------------------------------------------------------- # @@ -223,12 +224,13 @@ # significance is computed according to the Li & Ma expression for ON and # OFF Poisson measurements, see # `here `__. -# Since astropy convolution kernels only accept integers, we first convert -# our required size in degrees to int depending on our pixel size. # +# +# Also, since the off counts are obtained with a Ring banckground estimation, the we specify `correlate_off=False` +# to avoid over correlation. # Using a convolution radius of 0.04 degrees -estimator = ExcessMapEstimator(0.04 * u.deg, selection_optional=[]) +estimator = ExcessMapEstimator(0.04 * u.deg, selection_optional=[], correlate_off=False) lima_maps = estimator.run(stacked_on_off) significance_map = lima_maps["sqrt_ts"] @@ -238,28 +240,29 @@ fig, (ax1, ax2) = plt.subplots( figsize=(11, 4), subplot_kw={"projection": lima_maps.geom.wcs}, ncols=2 ) - ax1.set_title("Significance map") significance_map.plot(ax=ax1, add_cbar=True) - ax2.set_title("Excess map") excess_map.plot(ax=ax2, add_cbar=True) plt.show() - -###################################################################### +# ##################################################################### # It is often important to look at the significance distribution outside # the exclusion region to check that the background estimation is not # contaminated by gamma-ray events. This can be the case when exclusion # regions are not large enough. Typically, we expect the off distribution -# to be a standard normal distribution. +# to be a standard normal distribution. To compute the significance distribution outside the exclusion region, +# we can recompute the maps after adding a `mask_fit` to out dataset. # -# create a 2D mask for the images -significance_map_off = significance_map * exclusion_mask +# Mask the regions with gamma-ray emission +stacked_on_off.mask_fit = exclusion_mask +lima_maps2 = estimator.run(stacked_on_off) +significance_map_off = lima_maps2["sqrt_ts"] -kwargs_axes = {"xlabel": "Significance", "yscale": "log", "ylim": (1e-5, 1)} +# region +kwargs_axes = {"xlabel": "Significance", "yscale": "log", "ylim": (1e-3, 1)} ax, _ = plot_distribution( significance_map, kwargs_hist={ @@ -267,7 +270,7 @@ "alpha": 0.5, "color": "red", "label": "all bins", - "bins": 21, + "bins": 51, }, kwargs_axes=kwargs_axes, ) @@ -281,11 +284,12 @@ "alpha": 0.5, "color": "blue", "label": "off bins", - "bins": 21, + "bins": 51, }, kwargs_axes=kwargs_axes, ) plt.show() +# endregion # sphinx_gallery_thumbnail_number = 2 From dc79c23ff34c1926670bdaae342559d481088cdb Mon Sep 17 00:00:00 2001 From: Atreyee Sinha Date: Mon, 21 Oct 2024 18:13:49 +0200 Subject: [PATCH 19/55] Update examples/tutorials/analysis-2d/ring_background.py Co-authored-by: Daniel Morcuende --- examples/tutorials/analysis-2d/ring_background.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/tutorials/analysis-2d/ring_background.py b/examples/tutorials/analysis-2d/ring_background.py index 7ffb20e6df..7ec4dcea9d 100644 --- a/examples/tutorials/analysis-2d/ring_background.py +++ b/examples/tutorials/analysis-2d/ring_background.py @@ -226,7 +226,7 @@ # `here `__. # # -# Also, since the off counts are obtained with a Ring banckground estimation, the we specify `correlate_off=False` +# Also, since the off counts are obtained with a ring background estimation, then we specify `correlate_off=False` # to avoid over correlation. # Using a convolution radius of 0.04 degrees From 778140189afe6bbd9e3e833f309bd9922f8c364b Mon Sep 17 00:00:00 2001 From: Atreyee Sinha Date: Mon, 21 Oct 2024 18:13:58 +0200 Subject: [PATCH 20/55] Update examples/tutorials/analysis-2d/ring_background.py Co-authored-by: Daniel Morcuende --- examples/tutorials/analysis-2d/ring_background.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/tutorials/analysis-2d/ring_background.py b/examples/tutorials/analysis-2d/ring_background.py index 7ec4dcea9d..50c5898643 100644 --- a/examples/tutorials/analysis-2d/ring_background.py +++ b/examples/tutorials/analysis-2d/ring_background.py @@ -252,7 +252,7 @@ # contaminated by gamma-ray events. This can be the case when exclusion # regions are not large enough. Typically, we expect the off distribution # to be a standard normal distribution. To compute the significance distribution outside the exclusion region, -# we can recompute the maps after adding a `mask_fit` to out dataset. +# we can recompute the maps after adding a `mask_fit` to our dataset. # # Mask the regions with gamma-ray emission From baa73f3451da682b39e863efd965258703e83055 Mon Sep 17 00:00:00 2001 From: Atreyee Sinha Date: Mon, 21 Oct 2024 18:17:09 +0200 Subject: [PATCH 21/55] Update examples/tutorials/analysis-2d/ring_background.py --- examples/tutorials/analysis-2d/ring_background.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/tutorials/analysis-2d/ring_background.py b/examples/tutorials/analysis-2d/ring_background.py index 50c5898643..f30db9a7db 100644 --- a/examples/tutorials/analysis-2d/ring_background.py +++ b/examples/tutorials/analysis-2d/ring_background.py @@ -39,7 +39,7 @@ """ -# ##################################################################### +##################################################################### # Setup # ----- # From d616eb3539eb1221417e87c41459cf958d4eb457 Mon Sep 17 00:00:00 2001 From: Atreyee Sinha Date: Mon, 21 Oct 2024 18:17:16 +0200 Subject: [PATCH 22/55] Update examples/tutorials/analysis-2d/ring_background.py --- examples/tutorials/analysis-2d/ring_background.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/tutorials/analysis-2d/ring_background.py b/examples/tutorials/analysis-2d/ring_background.py index f30db9a7db..09d01213ae 100644 --- a/examples/tutorials/analysis-2d/ring_background.py +++ b/examples/tutorials/analysis-2d/ring_background.py @@ -71,7 +71,7 @@ check_tutorials_setup() -# ##################################################################### +##################################################################### # Creating the config file # ------------------------ # From 47877982bdf4593c2fbbad9d179f5c55f928fb0f Mon Sep 17 00:00:00 2001 From: Atreyee Sinha Date: Mon, 21 Oct 2024 18:18:23 +0200 Subject: [PATCH 23/55] Update examples/tutorials/analysis-2d/ring_background.py --- examples/tutorials/analysis-2d/ring_background.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/tutorials/analysis-2d/ring_background.py b/examples/tutorials/analysis-2d/ring_background.py index 09d01213ae..ded57f7d7d 100644 --- a/examples/tutorials/analysis-2d/ring_background.py +++ b/examples/tutorials/analysis-2d/ring_background.py @@ -63,7 +63,7 @@ log = logging.getLogger(__name__) -# ##################################################################### +##################################################################### # Check setup # ----------- From a75fd4f8a85e8694e33801f8fbfea2cb7f9aeda2 Mon Sep 17 00:00:00 2001 From: Atreyee Sinha Date: Mon, 21 Oct 2024 18:18:32 +0200 Subject: [PATCH 24/55] Update examples/tutorials/analysis-2d/ring_background.py --- examples/tutorials/analysis-2d/ring_background.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/tutorials/analysis-2d/ring_background.py b/examples/tutorials/analysis-2d/ring_background.py index ded57f7d7d..885cdd2875 100644 --- a/examples/tutorials/analysis-2d/ring_background.py +++ b/examples/tutorials/analysis-2d/ring_background.py @@ -186,7 +186,7 @@ ) -# ##################################################################### +##################################################################### # Create a stacked dataset # ~~~~~~~~~~~~~~~~~~~~~~~~ # From 0557e174dcde26e2f96c1f4799f42c334053fb6d Mon Sep 17 00:00:00 2001 From: Atreyee Sinha Date: Mon, 21 Oct 2024 18:18:38 +0200 Subject: [PATCH 25/55] Update examples/tutorials/analysis-2d/ring_background.py --- examples/tutorials/analysis-2d/ring_background.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/tutorials/analysis-2d/ring_background.py b/examples/tutorials/analysis-2d/ring_background.py index 885cdd2875..daad272fe1 100644 --- a/examples/tutorials/analysis-2d/ring_background.py +++ b/examples/tutorials/analysis-2d/ring_background.py @@ -118,7 +118,7 @@ print(config) -# ##################################################################### +##################################################################### # Getting the reduced dataset # --------------------------- # From 4153630279f1d0564ef1369a26fb428656320262 Mon Sep 17 00:00:00 2001 From: Atreyee Sinha Date: Mon, 21 Oct 2024 18:18:45 +0200 Subject: [PATCH 26/55] Update examples/tutorials/analysis-2d/ring_background.py --- examples/tutorials/analysis-2d/ring_background.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/tutorials/analysis-2d/ring_background.py b/examples/tutorials/analysis-2d/ring_background.py index daad272fe1..1c057b6c87 100644 --- a/examples/tutorials/analysis-2d/ring_background.py +++ b/examples/tutorials/analysis-2d/ring_background.py @@ -155,7 +155,7 @@ # -# ##################################################################### +##################################################################### # Create exclusion mask # ~~~~~~~~~~~~~~~~~~~~~ # From 1085a9674cc1bc3db6869ecac5bc2c1578aff4d4 Mon Sep 17 00:00:00 2001 From: Atreyee Sinha Date: Mon, 21 Oct 2024 18:18:51 +0200 Subject: [PATCH 27/55] Update examples/tutorials/analysis-2d/ring_background.py --- examples/tutorials/analysis-2d/ring_background.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/tutorials/analysis-2d/ring_background.py b/examples/tutorials/analysis-2d/ring_background.py index 1c057b6c87..2490038593 100644 --- a/examples/tutorials/analysis-2d/ring_background.py +++ b/examples/tutorials/analysis-2d/ring_background.py @@ -176,7 +176,7 @@ plt.show() -# ##################################################################### +##################################################################### # For the present analysis, we use a ring with an inner radius of 0.5 deg # and width of 0.3 deg. # From 0947c2e928783f832a58471252597cd38186ebe1 Mon Sep 17 00:00:00 2001 From: Atreyee Sinha Date: Mon, 21 Oct 2024 18:18:57 +0200 Subject: [PATCH 28/55] Update examples/tutorials/analysis-2d/ring_background.py --- examples/tutorials/analysis-2d/ring_background.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/tutorials/analysis-2d/ring_background.py b/examples/tutorials/analysis-2d/ring_background.py index 2490038593..8e6049c32d 100644 --- a/examples/tutorials/analysis-2d/ring_background.py +++ b/examples/tutorials/analysis-2d/ring_background.py @@ -206,7 +206,7 @@ stacked_on_off.stack(dataset_on_off) -# ##################################################################### +##################################################################### # This `stacked_on_off` has `on` and `off` counts and acceptance # maps which we will use in all further analysis. The `acceptance` and # `acceptance_off` maps are the system acceptance of gamma-ray like From e979ab97cbabf85ccee20fd1ce44436bd1969dc8 Mon Sep 17 00:00:00 2001 From: Atreyee Sinha Date: Mon, 21 Oct 2024 18:19:03 +0200 Subject: [PATCH 29/55] Update examples/tutorials/analysis-2d/ring_background.py --- examples/tutorials/analysis-2d/ring_background.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/tutorials/analysis-2d/ring_background.py b/examples/tutorials/analysis-2d/ring_background.py index 8e6049c32d..7a24c906f8 100644 --- a/examples/tutorials/analysis-2d/ring_background.py +++ b/examples/tutorials/analysis-2d/ring_background.py @@ -226,7 +226,7 @@ # `here `__. # # -# Also, since the off counts are obtained with a ring background estimation, then we specify `correlate_off=False` +# Also, since the off counts are obtained with a ring background estimation, and are thus already correlated, we specify `correlate_off=False` # to avoid over correlation. # Using a convolution radius of 0.04 degrees From e26eee1d0f0854e819a79ab8dcbd5c23e078b9e4 Mon Sep 17 00:00:00 2001 From: Atreyee Sinha Date: Mon, 21 Oct 2024 18:20:22 +0200 Subject: [PATCH 30/55] Update examples/tutorials/analysis-2d/ring_background.py --- examples/tutorials/analysis-2d/ring_background.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/tutorials/analysis-2d/ring_background.py b/examples/tutorials/analysis-2d/ring_background.py index 7a24c906f8..b640aba9b4 100644 --- a/examples/tutorials/analysis-2d/ring_background.py +++ b/examples/tutorials/analysis-2d/ring_background.py @@ -246,7 +246,7 @@ excess_map.plot(ax=ax2, add_cbar=True) plt.show() -# ##################################################################### +##################################################################### # It is often important to look at the significance distribution outside # the exclusion region to check that the background estimation is not # contaminated by gamma-ray events. This can be the case when exclusion From e1de9e65689ef86a44931d9614fbd2fba3fe3641 Mon Sep 17 00:00:00 2001 From: Astro-Kirsty Date: Mon, 21 Oct 2024 19:07:41 +0200 Subject: [PATCH 31/55] revert changes to titles Signed-off-by: Astro-Kirsty --- .../tutorials/analysis-2d/ring_background.py | 22 +++++++++---------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/examples/tutorials/analysis-2d/ring_background.py b/examples/tutorials/analysis-2d/ring_background.py index b640aba9b4..aa8a0642aa 100644 --- a/examples/tutorials/analysis-2d/ring_background.py +++ b/examples/tutorials/analysis-2d/ring_background.py @@ -39,7 +39,7 @@ """ -##################################################################### +###################################################################### # Setup # ----- # @@ -63,7 +63,7 @@ log = logging.getLogger(__name__) -##################################################################### +###################################################################### # Check setup # ----------- @@ -71,7 +71,7 @@ check_tutorials_setup() -##################################################################### +###################################################################### # Creating the config file # ------------------------ # @@ -118,7 +118,7 @@ print(config) -##################################################################### +###################################################################### # Getting the reduced dataset # --------------------------- # @@ -145,7 +145,7 @@ analysis.get_datasets() -# ##################################################################### +###################################################################### # Extracting the ring background # ------------------------------ # @@ -155,7 +155,7 @@ # -##################################################################### +###################################################################### # Create exclusion mask # ~~~~~~~~~~~~~~~~~~~~~ # @@ -176,7 +176,7 @@ plt.show() -##################################################################### +###################################################################### # For the present analysis, we use a ring with an inner radius of 0.5 deg # and width of 0.3 deg. # @@ -186,7 +186,7 @@ ) -##################################################################### +###################################################################### # Create a stacked dataset # ~~~~~~~~~~~~~~~~~~~~~~~~ # @@ -206,7 +206,7 @@ stacked_on_off.stack(dataset_on_off) -##################################################################### +###################################################################### # This `stacked_on_off` has `on` and `off` counts and acceptance # maps which we will use in all further analysis. The `acceptance` and # `acceptance_off` maps are the system acceptance of gamma-ray like @@ -216,7 +216,7 @@ print(stacked_on_off) -# ##################################################################### +###################################################################### # Compute correlated significance and correlated excess maps # ---------------------------------------------------------- # @@ -246,7 +246,7 @@ excess_map.plot(ax=ax2, add_cbar=True) plt.show() -##################################################################### +###################################################################### # It is often important to look at the significance distribution outside # the exclusion region to check that the background estimation is not # contaminated by gamma-ray events. This can be the case when exclusion From 32ea414be94eba82368783bb2c58b1b1b99d4b2e Mon Sep 17 00:00:00 2001 From: Daniel Morcuende Date: Mon, 21 Oct 2024 13:18:03 +0200 Subject: [PATCH 32/55] Use float in PSF tick formatter to unify it with the one used in PSFMap --- gammapy/irf/psf/core.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/gammapy/irf/psf/core.py b/gammapy/irf/psf/core.py index 226d4ea594..e55cd7ffd4 100644 --- a/gammapy/irf/psf/core.py +++ b/gammapy/irf/psf/core.py @@ -168,7 +168,7 @@ def plot_containment_radius_vs_energy( ax.set_ylabel( f"Containment radius [{ax.yaxis.units.to_string(UNIT_STRING_FORMAT)}]" ) - ax.yaxis.set_major_formatter(mtick.FormatStrFormatter("%.1e")) + ax.yaxis.set_major_formatter(mtick.FormatStrFormatter("%.2f")) return ax def plot_containment_radius( From 6a6de469c245e02df39995d1e50399b91b384ded Mon Sep 17 00:00:00 2001 From: Atreyee Date: Tue, 22 Oct 2024 15:31:12 +0200 Subject: [PATCH 33/55] implement comments Signed-off-by: Atreyee Signed-off-by: --- .../analysis-time/light_curve_flare.py | 18 +++++++++++------- 1 file changed, 11 insertions(+), 7 deletions(-) diff --git a/examples/tutorials/analysis-time/light_curve_flare.py b/examples/tutorials/analysis-time/light_curve_flare.py index 679750c25e..2f2e4ffa28 100644 --- a/examples/tutorials/analysis-time/light_curve_flare.py +++ b/examples/tutorials/analysis-time/light_curve_flare.py @@ -127,7 +127,7 @@ # Define time intervals # --------------------- # -# We create the list of time intervals, each of duration 10 mins. Each time interval is an +# We create the list of time intervals, each of duration 10 minutes. Each time interval is an # `astropy.time.Time` object, containing a start and stop time. t0 = Time("2006-07-29T20:30") @@ -237,15 +237,16 @@ ###################################################################### # Define underlying model -# ---------------- +# ----------------------- # # Since we use forward folding to obtain the flux points in each bin, exact values will depend on the underlying model. In this example, we use a power law as used in the # `reference -# paper `. +# paper `__. # -# As we have are only using spectral datasats, we do not need any spatial models. +# As we have are only using spectral datasets, we do not need any spatial models. # -# **Note** : All time bins must have the same spectral model. To see how to investigate spectral variability, see :doc:`time resolved spectroscopy notebook `. +# **Note** : All time bins must have the same spectral model. To see how to investigate spectral variability, +# see :doc:`time resolved spectroscopy notebook `. spectral_model = PowerLawSpectralModel(amplitude=1e-10 * u.Unit("1 / (cm2 s TeV)")) sky_model = SkyModel(spatial_model=None, spectral_model=spectral_model, name="pks2155") @@ -271,7 +272,8 @@ # # We first create the `~gammapy.estimators.LightCurveEstimator` for the # list of datasets we just produced. We give the estimator the name of the -# source component to be fitted. We can directly compute the light curve in multiple energy bins by supplying a list of `energy_edges`. +# source component to be fitted. We can directly compute the light curve in multiple energy +# bins by supplying a list of `energy_edges`. # # By default the likelihood scan is computed from 0.2 to 5.0. # Here, we increase the max value to 10.0, because we are @@ -329,4 +331,6 @@ plt.legend() plt.show() -# We can use the sliced lightcurve to understand the variability, as shown in the doc:`/tutorials/analysis-time/variability_estimation` tutorial. +##################################################################### +# We can use the sliced lightcurve to understand the variability, +# as shown in the doc:`/tutorials/analysis-time/variability_estimation` tutorial. From ecacfdb27e706c7d08872753cc84d10b9f0f756c Mon Sep 17 00:00:00 2001 From: Atreyee Date: Tue, 22 Oct 2024 15:55:05 +0200 Subject: [PATCH 34/55] modify truncated x-axis Signed-off-by: Atreyee Signed-off-by: --- examples/tutorials/analysis-time/light_curve_flare.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/examples/tutorials/analysis-time/light_curve_flare.py b/examples/tutorials/analysis-time/light_curve_flare.py index 2f2e4ffa28..937dc5584e 100644 --- a/examples/tutorials/analysis-time/light_curve_flare.py +++ b/examples/tutorials/analysis-time/light_curve_flare.py @@ -304,6 +304,8 @@ # Finally we plot the result for the 1D lightcurve: # plt.figure(figsize=(8, 6)) +plt.tight_layout() +plt.subplots_adjust(bottom=0.3) lc_1d.plot(marker="o", axis_name="time", sed_type="flux") plt.show() @@ -326,6 +328,8 @@ lc_new = sliced_lc.resample_axis(axis_new) plt.figure(figsize=(8, 6)) +plt.tight_layout() +plt.subplots_adjust(bottom=0.3) ax = sliced_lc.plot(label="original") lc_new.plot(ax=ax, label="rebinned") plt.legend() From a6ce4be691bcc5242acec8840617034f5a592e1d Mon Sep 17 00:00:00 2001 From: Atreyee Date: Tue, 22 Oct 2024 16:22:55 +0200 Subject: [PATCH 35/55] add docs Signed-off-by: Atreyee Signed-off-by: --- docs/user-guide/howto.rst | 11 +++++++++++ examples/tutorials/api/estimators.py | 4 ++-- 2 files changed, 13 insertions(+), 2 deletions(-) diff --git a/docs/user-guide/howto.rst b/docs/user-guide/howto.rst index a058d500af..12a7d08abf 100644 --- a/docs/user-guide/howto.rst +++ b/docs/user-guide/howto.rst @@ -278,6 +278,17 @@ warning like so: .. accordion-footer:: +.. accordion-header:: + :id: HowToAvoidNaninFp + :title: Avoid `nan` results in Flux Point estimation + :link: ../tutorials/api/estimators.html#a-fully-configured-flux-points-estimation + +Sometimes, upper limit values may show as `nan` while running a `FluxPointsEstimator` or a `LightCurveEstimator`. +This often arises because the range of the norm parameter being scanned over is not sufficient. Increasing this range +usually solves the problem. In some cases, you can also consider configuring the estimator with a different `Fit` backend. + +.. accordion-footer:: + .. accordion-header:: :id: HowToProgressBar :title: Display a progress bar diff --git a/examples/tutorials/api/estimators.py b/examples/tutorials/api/estimators.py index 01cb717955..db7624f948 100644 --- a/examples/tutorials/api/estimators.py +++ b/examples/tutorials/api/estimators.py @@ -116,10 +116,10 @@ # fp_estimator.norm.scan_min = 0.1 -fp_estimator.norm.scan_max = 3 +fp_estimator.norm.scan_max = 10 ###################################################################### -# Note: The default scan range of the norm parameter is between 0.1 to 10. In case the upper +# Note: The default scan range of the norm parameter is between 0.1 to 5. In case the upper # limit values lie outside this range, nan values will be returned. It may thus be useful to # increase this range, specially for the computation of upper limits from weak sources. # From ccb65043ab77e0169804caf917ba14999096a9d5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Bruno=20Kh=C3=A9lifi?= Date: Tue, 22 Oct 2024 17:29:49 +0200 Subject: [PATCH 36/55] Update docs/user-guide/howto.rst Co-authored-by: Kirsty Feijen --- docs/user-guide/howto.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/user-guide/howto.rst b/docs/user-guide/howto.rst index 12a7d08abf..b23ee39738 100644 --- a/docs/user-guide/howto.rst +++ b/docs/user-guide/howto.rst @@ -280,7 +280,7 @@ warning like so: .. accordion-header:: :id: HowToAvoidNaninFp - :title: Avoid `nan` results in Flux Point estimation + :title: Avoid NaN results in Flux Point estimation :link: ../tutorials/api/estimators.html#a-fully-configured-flux-points-estimation Sometimes, upper limit values may show as `nan` while running a `FluxPointsEstimator` or a `LightCurveEstimator`. From 693930f93fc8f5599b202703f2488c268ba34f98 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Bruno=20Kh=C3=A9lifi?= Date: Tue, 22 Oct 2024 17:32:40 +0200 Subject: [PATCH 37/55] Update examples/tutorials/api/estimators.py Co-authored-by: Kirsty Feijen --- examples/tutorials/api/estimators.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/tutorials/api/estimators.py b/examples/tutorials/api/estimators.py index db7624f948..2255c165a5 100644 --- a/examples/tutorials/api/estimators.py +++ b/examples/tutorials/api/estimators.py @@ -119,7 +119,7 @@ fp_estimator.norm.scan_max = 10 ###################################################################### -# Note: The default scan range of the norm parameter is between 0.1 to 5. In case the upper +# Note: The default scan range of the norm parameter is between 0.2 to 5. In case the upper # limit values lie outside this range, nan values will be returned. It may thus be useful to # increase this range, specially for the computation of upper limits from weak sources. # From f95320208acfc5dd9d559934dc05db59a1b07082 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Bruno=20Kh=C3=A9lifi?= Date: Wed, 23 Oct 2024 10:33:50 +0200 Subject: [PATCH 38/55] Update examples/tutorials/analysis-time/light_curve_flare.py Co-authored-by: Kirsty Feijen --- examples/tutorials/analysis-time/light_curve_flare.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/tutorials/analysis-time/light_curve_flare.py b/examples/tutorials/analysis-time/light_curve_flare.py index 937dc5584e..315f856df9 100644 --- a/examples/tutorials/analysis-time/light_curve_flare.py +++ b/examples/tutorials/analysis-time/light_curve_flare.py @@ -312,7 +312,7 @@ ###################################################################### # Light curves once obtained can be rebinned using the likelihood profiles. -# Here, we rebin 3 adjacent bins together, to get 30 min bins. +# Here, we rebin 3 adjacent bins together, to get 30 minute bins. # # We will first slice `lc_1d` to obtain the lightcurve in the first energy bin # From e4e913f00bc17f71cebc4c20833e2716404b1908 Mon Sep 17 00:00:00 2001 From: Maxime Regeard Date: Wed, 23 Oct 2024 17:54:52 +0200 Subject: [PATCH 39/55] nomalise_phasecurve --- gammapy/modeling/models/temporal.py | 12 ++++++++++-- gammapy/modeling/models/tests/test_temporal.py | 14 ++++++-------- 2 files changed, 16 insertions(+), 10 deletions(-) diff --git a/gammapy/modeling/models/temporal.py b/gammapy/modeling/models/temporal.py index 39988f9dab..9bb498a72a 100644 --- a/gammapy/modeling/models/temporal.py +++ b/gammapy/modeling/models/temporal.py @@ -1,5 +1,6 @@ # Licensed under a 3-clause BSD style license - see LICENSE.rst """Time-dependent models.""" + import logging import numpy as np import scipy.interpolate @@ -7,7 +8,6 @@ from astropy.io import fits from astropy.table import Table from astropy.time import Time -from astropy.utils import lazyproperty from gammapy.maps import MapAxis, RegionNDMap, TimeMapAxis from gammapy.modeling import Parameter from gammapy.utils.compat import COPY_IF_NEEDED @@ -988,6 +988,14 @@ def __init__(self, table, filename=None, **kwargs): filename = str(make_path(filename)) self.filename = filename super().__init__(**kwargs) + self._normalise() + + def _normalise(self): + time_ref = Time(self.t_ref.value, scale=self.scale, format="mjd") + + norm_factor = 1 / (self.integral(time_ref, time_ref + 1 / self.f0.quantity)) + + self.table["NORM"] = norm_factor * self.table["NORM"] @classmethod def read( @@ -1065,7 +1073,7 @@ def write(self, path=None, overwrite=False): self.filename = str(make_path(path)) self.table.write(self.filename, overwrite=overwrite) - @lazyproperty + @property def _interpolator(self): x = self.table["PHASE"].data y = self.table["NORM"].data diff --git a/gammapy/modeling/models/tests/test_temporal.py b/gammapy/modeling/models/tests/test_temporal.py index 7ce68ab444..d4055609df 100644 --- a/gammapy/modeling/models/tests/test_temporal.py +++ b/gammapy/modeling/models/tests/test_temporal.py @@ -354,7 +354,6 @@ def test_sine_temporal_model_integral(): @requires_data() def test_to_dict(light_curve): - out = light_curve.to_dict() assert out["temporal"]["type"] == "LightCurveTemplateTemporalModel" assert "lightcrv_PKSB1222+216.fits" in out["temporal"]["filename"] @@ -362,7 +361,6 @@ def test_to_dict(light_curve): @requires_data() def test_with_skymodel(light_curve): - sky_model = SkyModel(spectral_model=PowerLawSpectralModel()) out = sky_model.to_dict() assert "temporal" not in out @@ -399,7 +397,7 @@ def test_phase_curve_model(tmp_path): ) result = phase_model(t_ref + [0, 0.025, 0.05] * u.s) - assert_allclose(result, [0, 0.5, 0], atol=1e-5) + assert_allclose(result, [0.000000e00, 1.999989e00, 2.169609e-05], atol=1e-5) phase_model.filename = str(make_path(tmp_path / "tmp.fits")) phase_model.write() @@ -415,17 +413,17 @@ def test_phase_curve_model(tmp_path): ) assert_allclose(new_model.table["PHASE"].data, phase) - assert_allclose(new_model.table["NORM"].data, norm) + assert_allclose(new_model.table["NORM"].data, norm * 4) # exact number of phases integral = phase_model.integral(t_ref, t_ref + 10 * u.s) - assert_allclose(integral, 0.25, rtol=1e-5) + assert_allclose(integral, 1.0, rtol=1e-5) # long duration. Should be equal to the phase average integral = phase_model.integral(t_ref + 1 * u.h, t_ref + 3 * u.h) - assert_allclose(integral, 0.25, rtol=1e-5) + assert_allclose(integral, 1.0, rtol=1e-5) # 1.25 phase integral = phase_model.integral(t_ref, t_ref + 62.5 * u.ms) - assert_allclose(integral, 0.225, rtol=1e-5) + assert_allclose(integral, 1.0, rtol=1e-5) def test_phase_curve_model_sample_time(): @@ -472,7 +470,7 @@ def test_phasecurve_DC1(): times = Time(t_ref, format="mjd") + [0.0, 0.5, 0.65, 1.0] * P0 norm = model(times) - assert_allclose(norm, [0.05, 0.15, 1.0, 0.05]) + assert_allclose(norm, [0.294118, 0.882353, 5.882353, 0.294118], atol=1e-5) with mpl_plot_check(): model.plot_phasogram(n_points=200) From 82b3d90a9caafa3361ef512af6ecf67984e83be0 Mon Sep 17 00:00:00 2001 From: Maxime Regeard Date: Thu, 24 Oct 2024 15:23:35 +0200 Subject: [PATCH 40/55] normalise table column at init --- gammapy/modeling/models/temporal.py | 44 +++++++++++++++++++++-------- 1 file changed, 33 insertions(+), 11 deletions(-) diff --git a/gammapy/modeling/models/temporal.py b/gammapy/modeling/models/temporal.py index 9bb498a72a..46253a582a 100644 --- a/gammapy/modeling/models/temporal.py +++ b/gammapy/modeling/models/temporal.py @@ -9,6 +9,7 @@ from astropy.table import Table from astropy.time import Time from gammapy.maps import MapAxis, RegionNDMap, TimeMapAxis +from astropy.utils import lazyproperty from gammapy.modeling import Parameter from gammapy.utils.compat import COPY_IF_NEEDED from gammapy.utils.random import InverseCDFSampler, get_random_state @@ -983,19 +984,30 @@ class TemplatePhaseCurveTemporalModel(TemporalModel): f2 = Parameter("f2", _f2_default, frozen=True) def __init__(self, table, filename=None, **kwargs): - self.table = table + self.table = self._normalise_table(table) if filename is not None: filename = str(make_path(filename)) self.filename = filename super().__init__(**kwargs) - self._normalise() - def _normalise(self): - time_ref = Time(self.t_ref.value, scale=self.scale, format="mjd") + @staticmethod + def _normalise_table(table): + x = table["PHASE"].data + y = table["NORM"].data + + try: + w = table["WEIGHT"].data + except KeyError: + w = None - norm_factor = 1 / (self.integral(time_ref, time_ref + 1 / self.f0.quantity)) + interpolator = scipy.interpolate.InterpolatedUnivariateSpline( + x, y, w=w, k=1, ext=2, bbox=[0.0, 1.0] + ) - self.table["NORM"] = norm_factor * self.table["NORM"] + integral = interpolator.integral(0, 1) + + table["NORM"] *= 1 / integral + return table @classmethod def read( @@ -1018,6 +1030,7 @@ def read( Filename with path. """ filename = str(make_path(path)) + return cls( Table.read(filename), filename=filename, @@ -1073,13 +1086,18 @@ def write(self, path=None, overwrite=False): self.filename = str(make_path(path)) self.table.write(self.filename, overwrite=overwrite) - @property + @lazyproperty def _interpolator(self): x = self.table["PHASE"].data y = self.table["NORM"].data + try: + w = self.table["WEIGHT"].data + except KeyError: + w = None + return scipy.interpolate.InterpolatedUnivariateSpline( - x, y, k=1, ext=2, bbox=[0.0, 1.0] + x, y, w=w, k=1, ext=2, bbox=[0.0, 1.0] ) def evaluate(self, time, t_ref, phi_ref, f0, f1, f2): @@ -1126,11 +1144,15 @@ def integral(self, t_min, t_max): ) - self._interpolator.antiderivative()(ph_min) # Divide by Jacobian (here we neglect variations of frequency during the integration period) - total = (phase_integral + start_integral + end_integral) / frequency + total = phase_integral + start_integral + end_integral # Normalize by total integration time - integral_norm = total / self.time_sum(t_min, t_max) + n_period = (self.time_sum(t_min, t_max) * frequency).to("") + if int(n_period) == 0: + n_period = 1 + + integral_norm = total / n_period - return integral_norm.to("") + return integral_norm @classmethod def from_dict(cls, data): From b40549d474f0b1f7cead5e67b98ad5141d71d237 Mon Sep 17 00:00:00 2001 From: Maxime Regeard Date: Thu, 24 Oct 2024 15:24:00 +0200 Subject: [PATCH 41/55] adapt test --- gammapy/modeling/models/tests/test_temporal.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/gammapy/modeling/models/tests/test_temporal.py b/gammapy/modeling/models/tests/test_temporal.py index d4055609df..2a726f8cbd 100644 --- a/gammapy/modeling/models/tests/test_temporal.py +++ b/gammapy/modeling/models/tests/test_temporal.py @@ -423,7 +423,7 @@ def test_phase_curve_model(tmp_path): assert_allclose(integral, 1.0, rtol=1e-5) # 1.25 phase integral = phase_model.integral(t_ref, t_ref + 62.5 * u.ms) - assert_allclose(integral, 1.0, rtol=1e-5) + assert_allclose(integral, 0.9, rtol=1e-5) def test_phase_curve_model_sample_time(): From 2ed9346ead8068776c5a722b7d1a354a4b9ded11 Mon Sep 17 00:00:00 2001 From: Maxime Regeard Date: Thu, 24 Oct 2024 15:55:15 +0200 Subject: [PATCH 42/55] float test --- gammapy/modeling/models/tests/test_temporal.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/gammapy/modeling/models/tests/test_temporal.py b/gammapy/modeling/models/tests/test_temporal.py index 2a726f8cbd..c8e78676a3 100644 --- a/gammapy/modeling/models/tests/test_temporal.py +++ b/gammapy/modeling/models/tests/test_temporal.py @@ -428,7 +428,7 @@ def test_phase_curve_model(tmp_path): def test_phase_curve_model_sample_time(): phase = np.linspace(0.0, 1, 51) - norm = 1 * (phase < 0.5) + norm = 1.0 * (phase < 0.5) table = Table(data={"PHASE": phase, "NORM": norm}) t_ref = Time("2020-06-01", scale="utc") From 89ce0a1805340123ba9af3deae00ff810bd7e05b Mon Sep 17 00:00:00 2001 From: Astro-Kirsty Date: Thu, 24 Oct 2024 17:55:43 +0200 Subject: [PATCH 43/55] Fix HowTo accordions Signed-off-by: Astro-Kirsty --- gammapy/utils/docs.py | 29 +++++++++++++++-------------- 1 file changed, 15 insertions(+), 14 deletions(-) diff --git a/gammapy/utils/docs.py b/gammapy/utils/docs.py index 257119c7c9..9a27e7c1ce 100644 --- a/gammapy/utils/docs.py +++ b/gammapy/utils/docs.py @@ -12,10 +12,11 @@ - https://doughellmann.com/blog/2010/05/09/defining-custom-roles-in-sphinx/ - http://docutils.sourceforge.net/docs/howto/rst-directives.html -- https://github.com/docutils-mirror/docutils/blob/master/docutils/parsers/rst/directives/images.py +- https://github.com/docutils/docutils/blob/master/docutils/docutils/parsers/rst/directives/images.py - https://github.com/sphinx-doc/sphinx/blob/master/sphinx/directives/other.py -- https://github.com/bokeh/bokeh/tree/master/bokeh/sphinxext +- https://github.com/bokeh/bokeh/tree/master/src/bokeh/sphinxext """ + import os from pathlib import Path from docutils.parsers.rst import Directive @@ -44,23 +45,23 @@ def run(self): raw = f"""
-