.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "tutorials/data/hawc.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. or to run this example in your browser via Binder .. rst-class:: sphx-glr-example-title .. _sphx_glr_tutorials_data_hawc.py: HAWC with Gammapy ================= Explore HAWC event lists and instrument response functions (IRFs), then perform the data reduction steps. Introduction ------------ `HAWC `__ is an array of water Cherenkov detectors located in Mexico that detects gamma-rays in the range between hundreds of GeV and hundreds of TeV. Gammapy recently added support of HAWC high level data analysis, after export to the current `open data level 3 format `__. The HAWC data is largely private. However, in 2022, a small sub-set of HAWC Pass4 event lists from the Crab nebula region was publicly `released `__. This dataset is 1 GB in size, so only a subset will be used here as an example. Tutorial overview ----------------- This notebook is a quick introduction to HAWC data analysis with Gammapy. It briefly describes the HAWC data and how to access it, and then uses a subset of the data to produce a `~gammapy.datasets.MapDataset`, to show how the data reduction is performed. The HAWC data release contains events where the energy is estimated using two different algorithms, referred to as "NN" and "GP" (see this `paper `__ for a detailed description). These two event classes are not independent, meaning that events are repeated between the NN and GP datasets. Therefore, these data should never be analysed jointly, and one of the two estimators needs to be chosen before proceeding. Once the data has been reduced to a `~gammapy.datasets.MapDataset`, there are no differences in the way that HAWC data is handled with respect to data from any other observatory, such as H.E.S.S. or CTAO. HAWC data access and reduction ------------------------------ This is how to access data and IRFs from the HAWC Crab event data release. .. GENERATED FROM PYTHON SOURCE LINES 50-61 .. code-block:: Python import astropy.units as u from astropy.coordinates import SkyCoord import matplotlib.pyplot as plt from gammapy.data import DataStore, HDUIndexTable, ObservationTable from gammapy.datasets import MapDataset from gammapy.estimators import ExcessMapEstimator from gammapy.makers import MapDatasetMaker, SafeMaskMaker from gammapy.maps import Map, MapAxis, WcsGeom .. GENERATED FROM PYTHON SOURCE LINES 62-64 Check setup ----------- .. GENERATED FROM PYTHON SOURCE LINES 64-68 .. code-block:: Python from gammapy.utils.check import check_tutorials_setup check_tutorials_setup() .. rst-class:: sphx-glr-script-out .. code-block:: none System: python_executable : /home/runner/work/gammapy-docs/gammapy-docs/gammapy/.tox/build_docs/bin/python python_version : 3.9.20 machine : x86_64 system : Linux Gammapy package: version : 1.3 path : /home/runner/work/gammapy-docs/gammapy-docs/gammapy/.tox/build_docs/lib/python3.9/site-packages/gammapy Other packages: numpy : 1.26.4 scipy : 1.13.1 astropy : 5.2.2 regions : 0.8 click : 8.1.7 yaml : 6.0.2 IPython : 8.18.1 jupyterlab : not installed matplotlib : 3.9.2 pandas : not installed healpy : 1.17.3 iminuit : 2.30.1 sherpa : 4.16.1 naima : 0.10.0 emcee : 3.1.6 corner : 2.2.3 ray : 2.39.0 Gammapy environment variables: GAMMAPY_DATA : /home/runner/work/gammapy-docs/gammapy-docs/gammapy-datasets/1.3 .. GENERATED FROM PYTHON SOURCE LINES 69-70 Chose which estimator we will use .. GENERATED FROM PYTHON SOURCE LINES 70-74 .. code-block:: Python energy_estimator = "NN" .. GENERATED FROM PYTHON SOURCE LINES 75-90 A useful way to organize the relevant files are with index tables. The observation index table contains information on each particular observation, such as the run ID. The HDU index table has a row per relevant file (i.e., events, effective area, psf…) and contains the path to said file. The HAWC data is divided into different event types, classified using the fraction of the array that was triggered by an event, a quantity usually referred to as "fHit". These event types are fully independent, meaning that an event will have a unique event type identifier, which is usually a number indicating which fHit bin the event corresponds to. For this reason, a single HAWC observation has several HDU index tables associated to it, one per event type. In each table, there will be paths to a distinct event list file and IRFs. In the HAWC event data release, all the HDU tables are saved into the same FITS file, and can be accesses through the choice of the hdu index. .. GENERATED FROM PYTHON SOURCE LINES 93-95 Load the tables ~~~~~~~~~~~~~~~ .. GENERATED FROM PYTHON SOURCE LINES 95-100 .. code-block:: Python data_path = "$GAMMAPY_DATA/hawc/crab_events_pass4/" hdu_filename = f"hdu-index-table-{energy_estimator}-Crab.fits.gz" obs_filename = f"obs-index-table-{energy_estimator}-Crab.fits.gz" .. GENERATED FROM PYTHON SOURCE LINES 101-102 There is only one observation table .. GENERATED FROM PYTHON SOURCE LINES 102-104 .. code-block:: Python obs_table = ObservationTable.read(data_path + obs_filename) .. GENERATED FROM PYTHON SOURCE LINES 105-109 The remainder of this tutorial utilises just one fHit value, however, for a regular analysis you should combine all fHit bins. Here, we utilise fHit bin number 6. We start by reading the HDU index table of this fHit bin .. GENERATED FROM PYTHON SOURCE LINES 109-114 .. code-block:: Python fHit = 6 hdu_table = HDUIndexTable.read(data_path + hdu_filename, hdu=fHit) .. GENERATED FROM PYTHON SOURCE LINES 115-116 From the tables, we can create a `~gammapy.data.DataStore`. .. GENERATED FROM PYTHON SOURCE LINES 116-122 .. code-block:: Python data_store = DataStore(hdu_table=hdu_table, obs_table=obs_table) data_store.info() .. rst-class:: sphx-glr-script-out .. code-block:: none Data store: HDU index table: BASE_DIR: /home/runner/work/gammapy-docs/gammapy-docs/gammapy-datasets/1.3/hawc/crab_events_pass4 Rows: 6 OBS_ID: 103000133 -- 103000133 HDU_TYPE: ['aeff', 'bkg', 'edisp', 'events', 'gti', 'psf'] HDU_CLASS: ['edisp_kernel_map', 'events', 'gti', 'map', 'psf_map_reco'] Observation table: Observatory name: 'N/A' Number of observations: 1 .. GENERATED FROM PYTHON SOURCE LINES 123-124 There is only one observation .. GENERATED FROM PYTHON SOURCE LINES 124-127 .. code-block:: Python obs = data_store.get_observations()[0] .. GENERATED FROM PYTHON SOURCE LINES 128-129 Peek events from this observation .. GENERATED FROM PYTHON SOURCE LINES 129-134 .. code-block:: Python obs.events.peek() plt.show() .. image-sg:: /tutorials/data/images/sphx_glr_hawc_001.png :alt: hawc :srcset: /tutorials/data/images/sphx_glr_hawc_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 135-136 Peek the energy dispersion: .. GENERATED FROM PYTHON SOURCE LINES 136-140 .. code-block:: Python obs.edisp.peek() plt.show() .. image-sg:: /tutorials/data/images/sphx_glr_hawc_002.png :alt: hawc :srcset: /tutorials/data/images/sphx_glr_hawc_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 141-142 Peek the psf: .. GENERATED FROM PYTHON SOURCE LINES 142-145 .. code-block:: Python obs.psf.peek() plt.show() .. image-sg:: /tutorials/data/images/sphx_glr_hawc_003.png :alt: Exposure, Containment radius at 1 TeV, Containment radius at center of map, PSF at center of map :srcset: /tutorials/data/images/sphx_glr_hawc_003.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 146-147 Peek the background for one transit: .. GENERATED FROM PYTHON SOURCE LINES 147-151 .. code-block:: Python plt.figure() obs.bkg.reduce_over_axes().plot(add_cbar=True) plt.show() .. image-sg:: /tutorials/data/images/sphx_glr_hawc_004.png :alt: hawc :srcset: /tutorials/data/images/sphx_glr_hawc_004.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 152-153 Peek the effective exposure for one transit: .. GENERATED FROM PYTHON SOURCE LINES 153-158 .. code-block:: Python plt.figure() obs.aeff.reduce_over_axes().plot(add_cbar=True) plt.show() .. image-sg:: /tutorials/data/images/sphx_glr_hawc_005.png :alt: hawc :srcset: /tutorials/data/images/sphx_glr_hawc_005.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 159-165 Data reduction into a MapDataset -------------------------------- We will now produce a `~gammapy.datasets.MapDataset` using the data from one of the fHit bins. In order to use all bins, one just needs to repeat this process inside a for loop that modifies the variable fHit. .. GENERATED FROM PYTHON SOURCE LINES 168-169 First we define the geometry and axes, starting with the energy reco axis: .. GENERATED FROM PYTHON SOURCE LINES 169-176 .. code-block:: Python energy_axis = MapAxis.from_edges( [1.00, 1.78, 3.16, 5.62, 10.0, 17.8, 31.6, 56.2, 100, 177, 316] * u.TeV, name="energy", interp="log", ) .. GENERATED FROM PYTHON SOURCE LINES 177-180 Note: this axis is the one used to create the background model map. If different edges are used, the `~gammapy.makers.MapDatasetMaker` will interpolate between them, which might lead to unexpected behaviour. .. GENERATED FROM PYTHON SOURCE LINES 182-183 Define the energy true axis: .. GENERATED FROM PYTHON SOURCE LINES 183-188 .. code-block:: Python energy_axis_true = MapAxis.from_energy_bounds( 1e-3, 1e4, nbin=140, unit="TeV", name="energy_true" ) .. GENERATED FROM PYTHON SOURCE LINES 189-190 Finally, create a geometry around the Crab location: .. GENERATED FROM PYTHON SOURCE LINES 190-199 .. code-block:: Python geom = WcsGeom.create( skydir=SkyCoord(ra=83.63, dec=22.01, unit="deg", frame="icrs"), width=6 * u.deg, axes=[energy_axis], binsz=0.05, ) .. GENERATED FROM PYTHON SOURCE LINES 200-201 Define the makers we will use: .. GENERATED FROM PYTHON SOURCE LINES 201-206 .. code-block:: Python maker = MapDatasetMaker(selection=["counts", "background", "exposure", "edisp", "psf"]) safe_mask_maker = SafeMaskMaker(methods=["aeff-max"], aeff_percent=10) .. GENERATED FROM PYTHON SOURCE LINES 207-210 Create an empty `~gammapy.datasets.MapDataset`. The keyword ``reco_psf=True`` is needed because the HAWC PSF is derived in reconstructed energy. .. GENERATED FROM PYTHON SOURCE LINES 210-216 .. code-block:: Python dataset_empty = MapDataset.create( geom, energy_axis_true=energy_axis_true, name=f"fHit {fHit}", reco_psf=True ) dataset = maker.run(dataset_empty, obs) .. GENERATED FROM PYTHON SOURCE LINES 217-222 The livetime information is used by the `~gammapy.makers.SafeMaskMaker` to retrieve the effective area from the exposure. The HAWC effective area is computed for one source transit above 45º zenith, which is around 6h. Since the effective area condition used here is relative to the maximum, this value does not actually impact the result. .. GENERATED FROM PYTHON SOURCE LINES 222-227 .. code-block:: Python dataset.exposure.meta["livetime"] = "6 h" dataset = safe_mask_maker.run(dataset) .. GENERATED FROM PYTHON SOURCE LINES 228-233 Now we have a dataset that has background and exposure quantities for one single transit, but our dataset might comprise more. The number of transits can be derived using the good time intervals (GTI) stored with the event list. For convenience, the HAWC data release already included this quantity as a map. .. GENERATED FROM PYTHON SOURCE LINES 233-237 .. code-block:: Python transit_map = Map.read(data_path + "irfs/TransitsMap_Crab.fits.gz") transit_number = transit_map.get_by_coord(geom.center_skydir) .. GENERATED FROM PYTHON SOURCE LINES 238-239 Correct the background and exposure quantities: .. GENERATED FROM PYTHON SOURCE LINES 239-243 .. code-block:: Python dataset.background.data *= transit_number dataset.exposure.data *= transit_number .. GENERATED FROM PYTHON SOURCE LINES 244-249 Check the dataset we produced ----------------------------- We will now check the contents of the dataset. We can use the ``.peek()`` method to quickly get a glimpse of the contents .. GENERATED FROM PYTHON SOURCE LINES 249-253 .. code-block:: Python dataset.peek() plt.show() .. image-sg:: /tutorials/data/images/sphx_glr_hawc_006.png :alt: Counts, Excess counts, Exposure, Background :srcset: /tutorials/data/images/sphx_glr_hawc_006.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 254-255 Create significance maps to check that the Crab is visible: .. GENERATED FROM PYTHON SOURCE LINES 255-266 .. code-block:: Python excess_estimator = ExcessMapEstimator( correlation_radius="0.2 deg", selection_optional=[], energy_edges=energy_axis.edges ) excess = excess_estimator.run(dataset) (dataset.mask * excess["sqrt_ts"]).plot_grid( add_cbar=True, cmap="coolwarm", vmin=-5, vmax=5 ) plt.show() .. image-sg:: /tutorials/data/images/sphx_glr_hawc_007.png :alt: Energy 1.00 TeV - 1.78 TeV, Energy 1.78 TeV - 3.16 TeV, Energy 3.16 TeV - 5.62 TeV, Energy 5.62 TeV - 10.0 TeV, Energy 10.0 TeV - 17.8 TeV, Energy 17.8 TeV - 31.6 TeV, Energy 31.6 TeV - 56.2 TeV, Energy 56.2 TeV - 100 TeV, Energy 100 TeV - 177 TeV, Energy 177 TeV - 316 TeV :srcset: /tutorials/data/images/sphx_glr_hawc_007.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 267-268 Combining all energies .. GENERATED FROM PYTHON SOURCE LINES 268-278 .. code-block:: Python excess_estimator_integrated = ExcessMapEstimator( correlation_radius="0.2 deg", selection_optional=[] ) excess_integrated = excess_estimator_integrated.run(dataset) excess_integrated["sqrt_ts"].plot(add_cbar=True) plt.show() .. image-sg:: /tutorials/data/images/sphx_glr_hawc_008.png :alt: hawc :srcset: /tutorials/data/images/sphx_glr_hawc_008.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 279-286 Exercises --------- - Repeat the process for a different fHit bin. - Repeat the process for all the fHit bins provided in the data release and fit a model to the result. .. GENERATED FROM PYTHON SOURCE LINES 289-296 Next steps ---------- Now you know how to access and work with HAWC data. All other tutorials and documentation concerning 3D analysis techniques and the `~gammapy.datasets.MapDataset` object can be used from this step on. .. _sphx_glr_download_tutorials_data_hawc.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: binder-badge .. image:: images/binder_badge_logo.svg :target: https://mybinder.org/v2/gh/gammapy/gammapy-webpage/v1.3?urlpath=lab/tree/notebooks/1.3/tutorials/data/hawc.ipynb :alt: Launch binder :width: 150 px .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: hawc.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: hawc.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: hawc.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_