.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "tutorials/analysis-1d/extended_source_spectral_analysis.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_analysis-1d_extended_source_spectral_analysis.py: Spectral analysis of extended sources ===================================== Perform a spectral analysis of an extended source. Prerequisites ------------- - Understanding of spectral analysis techniques in classical Cherenkov astronomy. - Understanding the basic data reduction and modeling/fitting processes with the gammapy library API as shown in the tutorial :doc:`/tutorials/starting/analysis_2` Context ------- Many VHE sources in the Galaxy are extended. Studying them with a 1D spectral analysis is more complex than studying point sources. One often has to use complex (i.e. non-circular) regions and more importantly, one has to take into account the fact that the instrument response is non-uniform over the selected region. A typical example is given by the supernova remnant RX J1713-3935 which is nearly 1 degree in diameter. See the `following article `__. **Objective: Measure the spectrum of RX J1713-3945 in a 1 degree region fully enclosing it.** Proposed approach ----------------- We have seen in the general presentation of the spectrum extraction for point sources (see :doc:`/tutorials/analysis-1d/spectral_analysis` tutorial) that Gammapy uses specific datasets makers to first produce reduced spectral data and then to extract OFF measurements with reflected background techniques: the `~gammapy.makers.SpectrumDatasetMaker` and the `~gammapy.makers.ReflectedRegionsBackgroundMaker`. However, if the flag `use_region_center` is not set to `False`, the former simply computes the reduced IRFs at the center of the ON region (assumed to be circular). This is no longer valid for extended sources. To be able to compute average responses in the ON region, we can set `use_region_center=False` with the `~gammapy.makers.SpectrumDatasetMaker`, in which case the values of the IRFs are averaged over the entire region. In summary, we have to: - Define an ON region (a `~regions.SkyRegion`) fully enclosing the source we want to study. - Define a `~gammapy.maps.RegionGeom` with the ON region and the required energy range (in particular, beware of the true energy range). - Create the necessary makers : - the spectrum dataset maker : `~gammapy.makers.SpectrumDatasetMaker` with `use_region_center=False` - the OFF background maker, here a `~gammapy.makers.ReflectedRegionsBackgroundMaker` - and usually the safe range maker : `~gammapy.makers.SafeMaskMaker` - Perform the data reduction loop. And for every observation: - Produce a spectrum dataset - Extract the OFF data to produce a `~gammapy.datasets.SpectrumDatasetOnOff` and compute a safe range for it. - Stack or store the resulting spectrum dataset. - Finally proceed with model fitting on the dataset as usual. Here, we will use the RX J1713-3945 observations from the H.E.S.S. first public test data release. The tutorial is implemented with the intermediate level API. Setup ----- As usual, we’ll start with some general imports… .. GENERATED FROM PYTHON SOURCE LINES 86-105 .. code-block:: Python import astropy.units as u from astropy.coordinates import Angle, SkyCoord from regions import CircleSkyRegion # %matplotlib inline import matplotlib.pyplot as plt from IPython.display import display from gammapy.data import DataStore from gammapy.datasets import Datasets, SpectrumDataset from gammapy.makers import ( ReflectedRegionsBackgroundMaker, SafeMaskMaker, SpectrumDatasetMaker, ) from gammapy.maps import MapAxis, RegionGeom from gammapy.modeling import Fit from gammapy.modeling.models import PowerLawSpectralModel, SkyModel .. GENERATED FROM PYTHON SOURCE LINES 106-108 Check setup ----------- .. GENERATED FROM PYTHON SOURCE LINES 108-113 .. 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 : /Users/mregeard/Workspace/dev/code/gammapy/gammapy/.tox/build_docs/bin/python python_version : 3.11.10 machine : x86_64 system : Darwin Gammapy package: version : 1.3.dev1205+g00f44f94ac path : /Users/mregeard/Workspace/dev/code/gammapy/gammapy/.tox/build_docs/lib/python3.11/site-packages/gammapy Other packages: numpy : 1.26.4 scipy : 1.14.1 astropy : 5.2.2 regions : 0.10 click : 8.1.7 yaml : 6.0.2 IPython : 8.28.0 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.2 ray : 2.37.0 Gammapy environment variables: GAMMAPY_DATA : /Users/mregeard/Workspace/dev/code/gammapy/gammapy-data/ .. GENERATED FROM PYTHON SOURCE LINES 114-120 Select the data --------------- We first set the datastore and retrieve a few observations from our source. .. GENERATED FROM PYTHON SOURCE LINES 120-129 .. code-block:: Python datastore = DataStore.from_dir("$GAMMAPY_DATA/hess-dl3-dr1/") obs_ids = [20326, 20327, 20349, 20350, 20396, 20397] # In case you want to use all RX J1713 data in the HESS DR1 # other_ids=[20421, 20422, 20517, 20518, 20519, 20521, 20898, 20899, 20900] observations = datastore.get_observations(obs_ids) .. GENERATED FROM PYTHON SOURCE LINES 130-133 Prepare the datasets creation ----------------------------- .. GENERATED FROM PYTHON SOURCE LINES 136-143 Select the ON region ~~~~~~~~~~~~~~~~~~~~ Here we take a simple 1 degree circular region because it fits well with the morphology of RX J1713-3945. More complex regions could be used e.g. `~regions.EllipseSkyRegion` or `~regions.RectangleSkyRegion`. .. GENERATED FROM PYTHON SOURCE LINES 143-149 .. code-block:: Python target_position = SkyCoord(347.3, -0.5, unit="deg", frame="galactic") radius = Angle("0.5 deg") on_region = CircleSkyRegion(target_position, radius) .. GENERATED FROM PYTHON SOURCE LINES 150-161 Define the geometries ~~~~~~~~~~~~~~~~~~~~~ This part is especially important. - We have to define first energy axes. They define the axes of the resulting `~gammapy.datasets.SpectrumDatasetOnOff`. In particular, we have to be careful to the true energy axis: it has to cover a larger range than the reconstructed energy one. - Then we define the region geometry itself from the on region. .. GENERATED FROM PYTHON SOURCE LINES 161-173 .. code-block:: Python # The binning of the final spectrum is defined here. energy_axis = MapAxis.from_energy_bounds(0.1, 40.0, 10, unit="TeV") # Reduced IRFs are defined in true energy (i.e. not measured energy). energy_axis_true = MapAxis.from_energy_bounds( 0.05, 100, 30, unit="TeV", name="energy_true" ) geom = RegionGeom(on_region, axes=[energy_axis]) .. GENERATED FROM PYTHON SOURCE LINES 174-179 Create the makers ~~~~~~~~~~~~~~~~~ First we instantiate the target `~gammapy.datasets.SpectrumDataset`. .. GENERATED FROM PYTHON SOURCE LINES 179-186 .. code-block:: Python dataset_empty = SpectrumDataset.create( geom=geom, energy_axis_true=energy_axis_true, ) .. GENERATED FROM PYTHON SOURCE LINES 187-196 Now we create its associated maker. Here we need to produce, counts, exposure and edisp (energy dispersion) entries. PSF and IRF background are not needed, therefore we don’t compute them. **IMPORTANT**: Note that `use_region_center` is set to `False`. This is necessary so that the `~gammapy.makers.SpectrumDatasetMaker` considers the whole region in the IRF computation and not only the center. .. GENERATED FROM PYTHON SOURCE LINES 196-202 .. code-block:: Python maker = SpectrumDatasetMaker( selection=["counts", "exposure", "edisp"], use_region_center=False ) .. GENERATED FROM PYTHON SOURCE LINES 203-207 Now we create the OFF background maker for the spectra. If we have an exclusion region, we have to pass it here. We also define the safe range maker. .. GENERATED FROM PYTHON SOURCE LINES 207-212 .. code-block:: Python bkg_maker = ReflectedRegionsBackgroundMaker() safe_mask_maker = SafeMaskMaker(methods=["aeff-max"], aeff_percent=10) .. GENERATED FROM PYTHON SOURCE LINES 213-223 Perform the data reduction loop. -------------------------------- We can now run over selected observations. For each of them, we: - Create the `~gammapy.datasets.SpectrumDataset` - Compute the OFF via the reflected background method and create a `~gammapy.datasets.SpectrumDatasetOnOff` object - Run the safe mask maker on it - Add the `~gammapy.datasets.SpectrumDatasetOnOff` to the list. .. GENERATED FROM PYTHON SOURCE LINES 225-243 .. code-block:: Python datasets = Datasets() for obs in observations: # A SpectrumDataset is filled in this geometry dataset = maker.run(dataset_empty.copy(name=f"obs-{obs.obs_id}"), obs) # Define safe mask dataset = safe_mask_maker.run(dataset, obs) # Compute OFF dataset = bkg_maker.run(dataset, obs) # Append dataset to the list datasets.append(dataset) display(datasets.meta_table) .. rst-class:: sphx-glr-script-out .. code-block:: none NAME TYPE TELESCOP OBS_ID OBS_MODE RA_PNT ... MJDREFI MJDREFF TIMESYS GEOLON GEOLAT ALTITUDE deg ... --------- -------------------- -------- ------ -------- --------------- ... ------- ------- ------- ---------------- ----------------- ----------------- obs-20326 SpectrumDatasetOnOff HESS 20326 POINTING 259.29851667325 ... 0 0.0 tt 16.5002222222222 -23.2717777777778 1834.999999999783 obs-20327 SpectrumDatasetOnOff HESS 20327 POINTING 257.47731666009 ... 0 0.0 tt 16.5002222222222 -23.2717777777778 1834.999999999783 obs-20349 SpectrumDatasetOnOff HESS 20349 POINTING 259.29851667325 ... 0 0.0 tt 16.5002222222222 -23.2717777777778 1834.999999999783 obs-20350 SpectrumDatasetOnOff HESS 20350 POINTING 257.47731666009 ... 0 0.0 tt 16.5002222222222 -23.2717777777778 1834.999999999783 obs-20396 SpectrumDatasetOnOff HESS 20396 POINTING 258.38791666667 ... 0 0.0 tt 16.5002222222222 -23.2717777777778 1834.999999999783 obs-20397 SpectrumDatasetOnOff HESS 20397 POINTING 258.38791666667 ... 0 0.0 tt 16.5002222222222 -23.2717777777778 1834.999999999783 .. GENERATED FROM PYTHON SOURCE LINES 244-249 Explore the results ------------------- We can peek at the content of the spectrum datasets .. GENERATED FROM PYTHON SOURCE LINES 249-254 .. code-block:: Python datasets[0].peek() plt.show() .. image-sg:: /tutorials/analysis-1d/images/sphx_glr_extended_source_spectral_analysis_001.png :alt: Counts, Exposure, Energy Dispersion :srcset: /tutorials/analysis-1d/images/sphx_glr_extended_source_spectral_analysis_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 255-262 Cumulative excess and significance ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Finally, we can look at cumulative significance and number of excesses. This is done with the `info_table` method of `~gammapy.datasets.Datasets`. .. GENERATED FROM PYTHON SOURCE LINES 262-267 .. code-block:: Python info_table = datasets.info_table(cumulative=True) display(info_table) .. rst-class:: sphx-glr-script-out .. code-block:: none name counts excess sqrt_ts background ... stat_sum counts_off acceptance acceptance_off alpha ... ------- ------ --------------- ------------------ ---------------- ... ------------------ ---------- ---------- ------------------ ------------------- stacked 1216 170.5 4.159464335903991 1045.5 ... 43.26800161601427 2091 9.0 18.0 0.5 stacked 2339 270.5 4.7224642932828935 2068.5 ... 72.2340282463576 4137 18.0 36.0 0.5 stacked 3521 480.5 6.880790051412004 3040.5 ... 121.08402714166986 6081 27.0 54.0 0.5 stacked 4684 653.0 8.11478193157773 4031.0 ... 159.53811351626462 8062 36.0 72.0 0.5 stacked 5895 874.66650390625 9.869911175403269 5020.33349609375 ... 214.8627716890898 11030 45.0 98.86793518066406 0.45515263080596924 stacked 6985 993.166015625 10.251109146411467 5991.833984375 ... 238.19702399208626 12973 54.0 116.91612243652344 0.46186956763267517 .. GENERATED FROM PYTHON SOURCE LINES 268-269 And make the corresponding plots .. GENERATED FROM PYTHON SOURCE LINES 269-294 .. code-block:: Python fig, (ax_excess, ax_sqrt_ts) = plt.subplots(figsize=(10, 4), ncols=2, nrows=1) ax_excess.plot( info_table["livetime"].to("h"), info_table["excess"], marker="o", ls="none", ) ax_excess.set_title("Excess") ax_excess.set_xlabel("Livetime [h]") ax_excess.set_ylabel("Excess events") ax_sqrt_ts.plot( info_table["livetime"].to("h"), info_table["sqrt_ts"], marker="o", ls="none", ) ax_sqrt_ts.set_title("Sqrt(TS)") ax_sqrt_ts.set_xlabel("Livetime [h]") ax_sqrt_ts.set_ylabel("Sqrt(TS)") plt.show() .. image-sg:: /tutorials/analysis-1d/images/sphx_glr_extended_source_spectral_analysis_002.png :alt: Excess, Sqrt(TS) :srcset: /tutorials/analysis-1d/images/sphx_glr_extended_source_spectral_analysis_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 295-303 Perform spectral model fitting ------------------------------ Here we perform a joint fit. We first create the model, here a simple powerlaw, and assign it to every dataset in the `~gammapy.datasets.Datasets`. .. GENERATED FROM PYTHON SOURCE LINES 303-312 .. code-block:: Python spectral_model = PowerLawSpectralModel( index=2, amplitude=2e-11 * u.Unit("cm-2 s-1 TeV-1"), reference=1 * u.TeV ) model = SkyModel(spectral_model=spectral_model, name="RXJ 1713") datasets.models = [model] .. GENERATED FROM PYTHON SOURCE LINES 313-315 Now we can run the fit .. GENERATED FROM PYTHON SOURCE LINES 315-321 .. code-block:: Python fit_joint = Fit() result_joint = fit_joint.run(datasets=datasets) print(result_joint) .. rst-class:: sphx-glr-script-out .. code-block:: none OptimizeResult backend : minuit method : migrad success : True message : Optimization terminated successfully. nfev : 38 total stat : 52.79 CovarianceResult backend : minuit method : hesse success : True message : Hesse terminated successfully. .. GENERATED FROM PYTHON SOURCE LINES 322-327 Explore the fit results ~~~~~~~~~~~~~~~~~~~~~~~ First the fitted parameters values and their errors. .. GENERATED FROM PYTHON SOURCE LINES 327-331 .. code-block:: Python display(datasets.models.to_parameters_table()) .. rst-class:: sphx-glr-script-out .. code-block:: none model type name value unit error min max frozen link prior -------- ---- --------- ---------- -------------- --------- --- --- ------ ---- ----- RXJ 1713 index 2.1102e+00 6.129e-02 nan nan False RXJ 1713 amplitude 1.3576e-11 cm-2 s-1 TeV-1 9.757e-13 nan nan False RXJ 1713 reference 1.0000e+00 TeV 0.000e+00 nan nan True .. GENERATED FROM PYTHON SOURCE LINES 332-336 Then plot the fit result to compare measured and expected counts. Rather than plotting them for each individual dataset, we stack all datasets and plot the fit result on the result. .. GENERATED FROM PYTHON SOURCE LINES 336-346 .. code-block:: Python # First stack them all reduced = datasets.stack_reduce() # Assign the fitted model reduced.models = model # Plot the result ax_spectrum, ax_residuals = reduced.plot_fit() plt.show() .. image-sg:: /tutorials/analysis-1d/images/sphx_glr_extended_source_spectral_analysis_003.png :alt: extended source spectral analysis :srcset: /tutorials/analysis-1d/images/sphx_glr_extended_source_spectral_analysis_003.png :class: sphx-glr-single-img .. _sphx_glr_download_tutorials_analysis-1d_extended_source_spectral_analysis.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/main?urlpath=lab/tree/notebooks/dev/tutorials/analysis-1d/extended_source_spectral_analysis.ipynb :alt: Launch binder :width: 150 px .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: extended_source_spectral_analysis.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: extended_source_spectral_analysis.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: extended_source_spectral_analysis.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_