.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "tutorials/analysis-3d/cta_data_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-3d_cta_data_analysis.py: Basic image exploration and fitting =================================== Detect sources, produce a sky image and a spectrum using CTA-1DC data. Introduction ------------ **This notebook shows an example how to make a sky image and spectrum for simulated CTAO data with Gammapy.** The dataset we will use is three observation runs on the Galactic Center. This is a tiny (and thus quick to process and play with and learn) subset of the simulated CTAO dataset that was produced for the first data challenge in August 2017. .. GENERATED FROM PYTHON SOURCE LINES 21-26 Setup ----- As usual, we’ll start with some setup … .. GENERATED FROM PYTHON SOURCE LINES 26-58 .. code-block:: Python # Configure the logger, so that the spectral analysis # isn't so chatty about what it's doing. import logging import astropy.units as u from astropy.coordinates import SkyCoord from regions import CircleSkyRegion import matplotlib.pyplot as plt from IPython.display import display from gammapy.data import DataStore from gammapy.datasets import Datasets, FluxPointsDataset, MapDataset, SpectrumDataset from gammapy.estimators import FluxPointsEstimator, TSMapEstimator from gammapy.estimators.utils import find_peaks from gammapy.makers import ( MapDatasetMaker, ReflectedRegionsBackgroundMaker, SafeMaskMaker, SpectrumDatasetMaker, ) from gammapy.maps import MapAxis, RegionGeom, WcsGeom from gammapy.modeling import Fit from gammapy.modeling.models import ( GaussianSpatialModel, PowerLawSpectralModel, SkyModel, ) from gammapy.visualization import plot_npred_signal, plot_spectrum_datasets_off_regions from gammapy.utils.check import check_tutorials_setup logging.basicConfig() log = logging.getLogger("gammapy.spectrum") log.setLevel(logging.ERROR) .. GENERATED FROM PYTHON SOURCE LINES 59-61 Check setup ----------- .. GENERATED FROM PYTHON SOURCE LINES 61-66 .. code-block:: Python 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 67-76 Select observations ------------------- A Gammapy analysis usually starts by creating a `~gammapy.data.DataStore` and selecting observations. This is shown in detail in other notebooks (see e.g. the :doc:`/tutorials/starting/analysis_2` tutorial), here we choose three observations near the Galactic Center. .. GENERATED FROM PYTHON SOURCE LINES 76-96 .. code-block:: Python data_store = DataStore.from_dir("$GAMMAPY_DATA/cta-1dc/index/gps") # Just as a reminder: this is how to select observations # from astropy.coordinates import SkyCoord # table = data_store.obs_table # pos_obs = SkyCoord(table['GLON_PNT'], table['GLAT_PNT'], frame='galactic', unit='deg') # pos_target = SkyCoord(0, 0, frame='galactic', unit='deg') # offset = pos_target.separation(pos_obs).deg # mask = (1 < offset) & (offset < 2) # table = table[mask] # table.show_in_browser(jsviewer=True) obs_id = [110380, 111140, 111159] observations = data_store.get_observations(obs_id) obs_cols = ["OBS_ID", "GLON_PNT", "GLAT_PNT", "LIVETIME"] display(data_store.obs_table.select_obs_id(obs_id)[obs_cols]) .. rst-class:: sphx-glr-script-out .. code-block:: none OBS_ID GLON_PNT GLAT_PNT LIVETIME deg deg s ------ ------------------ ------------------ -------- 110380 359.9999912037958 -1.299995937905366 1764.0 111140 358.4999833830074 1.3000020211954284 1764.0 111159 1.5000056568267741 1.299940468335294 1764.0 .. GENERATED FROM PYTHON SOURCE LINES 97-106 Make sky images --------------- Define map geometry ~~~~~~~~~~~~~~~~~~~ Select the target position and define an ON region for the spectral analysis .. GENERATED FROM PYTHON SOURCE LINES 106-127 .. code-block:: Python axis = MapAxis.from_energy_bounds( 0.1, 10, nbin=10, unit="TeV", name="energy", ) axis_true = MapAxis.from_energy_bounds( 0.05, 20, nbin=20, name="energy_true", unit="TeV", ) geom = WcsGeom.create( skydir=(0, 0), npix=(500, 400), binsz=0.02, frame="galactic", axes=[axis] ) print(geom) .. rst-class:: sphx-glr-script-out .. code-block:: none WcsGeom axes : ['lon', 'lat', 'energy'] shape : (500, 400, 10) ndim : 3 frame : galactic projection : CAR center : 0.0 deg, 0.0 deg width : 10.0 deg x 8.0 deg wcs ref : 0.0 deg, 0.0 deg .. GENERATED FROM PYTHON SOURCE LINES 128-131 Compute images ~~~~~~~~~~~~~~ .. GENERATED FROM PYTHON SOURCE LINES 133-151 .. code-block:: Python stacked = MapDataset.create(geom=geom, energy_axis_true=axis_true) maker = MapDatasetMaker(selection=["counts", "background", "exposure", "psf"]) maker_safe_mask = SafeMaskMaker(methods=["offset-max"], offset_max=2.5 * u.deg) for obs in observations: cutout = stacked.cutout(obs.get_pointing_icrs(obs.tmid), width="5 deg") dataset = maker.run(cutout, obs) dataset = maker_safe_mask.run(dataset, obs) stacked.stack(dataset) # # The maps are cubes, with an energy axis. # Let's also make some images: # dataset_image = stacked.to_image() geom_image = dataset_image.geoms["geom"] .. rst-class:: sphx-glr-script-out .. code-block:: none /home/runner/work/gammapy-docs/gammapy-docs/gammapy/.tox/build_docs/lib/python3.9/site-packages/astropy/units/core.py:2097: UnitsWarning: '1/s/MeV/sr' did not parse as fits unit: Numeric factor not supported by FITS If this is meant to be a custom unit, define it with 'u.def_unit'. To have it recognized inside a file reader or other code, enable it with 'u.add_enabled_units'. For details, see https://docs.astropy.org/en/latest/units/combining_and_defining.html warnings.warn(msg, UnitsWarning) /home/runner/work/gammapy-docs/gammapy-docs/gammapy/.tox/build_docs/lib/python3.9/site-packages/astropy/units/core.py:2097: UnitsWarning: '1/s/MeV/sr' did not parse as fits unit: Numeric factor not supported by FITS If this is meant to be a custom unit, define it with 'u.def_unit'. To have it recognized inside a file reader or other code, enable it with 'u.add_enabled_units'. For details, see https://docs.astropy.org/en/latest/units/combining_and_defining.html warnings.warn(msg, UnitsWarning) /home/runner/work/gammapy-docs/gammapy-docs/gammapy/.tox/build_docs/lib/python3.9/site-packages/astropy/units/core.py:2097: UnitsWarning: '1/s/MeV/sr' did not parse as fits unit: Numeric factor not supported by FITS If this is meant to be a custom unit, define it with 'u.def_unit'. To have it recognized inside a file reader or other code, enable it with 'u.add_enabled_units'. For details, see https://docs.astropy.org/en/latest/units/combining_and_defining.html warnings.warn(msg, UnitsWarning) .. GENERATED FROM PYTHON SOURCE LINES 152-157 Show images ~~~~~~~~~~~ Let’s have a quick look at the images we computed … .. GENERATED FROM PYTHON SOURCE LINES 157-176 .. code-block:: Python fig, (ax1, ax2, ax3) = plt.subplots( figsize=(15, 5), ncols=3, subplot_kw={"projection": geom_image.wcs}, gridspec_kw={"left": 0.1, "right": 0.9}, ) ax1.set_title("Counts map") dataset_image.counts.smooth(2).plot(ax=ax1, vmax=5) ax2.set_title("Background map") dataset_image.background.plot(ax=ax2, vmax=5) ax3.set_title("Excess map") dataset_image.excess.smooth(3).plot(ax=ax3, vmax=2) plt.show() .. image-sg:: /tutorials/analysis-3d/images/sphx_glr_cta_data_analysis_001.png :alt: Counts map, Background map, Excess map :srcset: /tutorials/analysis-3d/images/sphx_glr_cta_data_analysis_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 177-184 Source Detection ---------------- Use the class `~gammapy.estimators.TSMapEstimator` and function `~gammapy.estimators.utils.find_peaks` to detect sources on the images. We search for 0.1 deg sigma gaussian sources in the dataset. .. GENERATED FROM PYTHON SOURCE LINES 184-198 .. code-block:: Python spatial_model = GaussianSpatialModel(sigma="0.05 deg") spectral_model = PowerLawSpectralModel(index=2) model = SkyModel(spatial_model=spatial_model, spectral_model=spectral_model) ts_image_estimator = TSMapEstimator( model, kernel_width="0.5 deg", selection_optional=[], downsampling_factor=2, sum_over_energy_groups=False, energy_edges=[0.1, 10] * u.TeV, ) .. GENERATED FROM PYTHON SOURCE LINES 199-208 .. code-block:: Python images_ts = ts_image_estimator.run(stacked) sources = find_peaks( images_ts["sqrt_ts"], threshold=5, min_distance="0.2 deg", ) display(sources) .. rst-class:: sphx-glr-script-out .. code-block:: none value x y ra dec deg deg ------ --- --- --------- --------- 36.125 252 197 266.42400 -29.00490 17.968 207 202 266.85900 -28.18386 13.13 186 199 267.16303 -27.85527 9.8672 373 205 264.79470 -30.97749 8.3583 298 169 266.42267 -30.08192 8.2689 308 187 265.94723 -30.06430 6.3229 90 209 268.07455 -26.10409 5.1489 87 226 267.78333 -25.87897 5.0163 239 167 267.16511 -29.09348 .. GENERATED FROM PYTHON SOURCE LINES 209-211 To get the position of the sources, simply .. GENERATED FROM PYTHON SOURCE LINES 211-214 .. code-block:: Python source_pos = SkyCoord(sources["ra"], sources["dec"]) print(source_pos) .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 215-217 Plot sources on top of significance sky image .. GENERATED FROM PYTHON SOURCE LINES 217-233 .. code-block:: Python fig, ax = plt.subplots(figsize=(8, 6), subplot_kw={"projection": geom_image.wcs}) images_ts["sqrt_ts"].plot(ax=ax, add_cbar=True) ax.scatter( source_pos.ra.deg, source_pos.dec.deg, transform=ax.get_transform("icrs"), color="none", edgecolor="white", marker="o", s=200, lw=1.5, ) plt.show() .. image-sg:: /tutorials/analysis-3d/images/sphx_glr_cta_data_analysis_002.png :alt: cta data analysis :srcset: /tutorials/analysis-3d/images/sphx_glr_cta_data_analysis_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 234-239 Spatial analysis ---------------- See other notebooks for how to run a 3D cube or 2D image based analysis. .. GENERATED FROM PYTHON SOURCE LINES 242-249 Spectrum -------- We’ll run a spectral analysis using the classical reflected regions background estimation method, and using the on-off (often called WSTAT) likelihood function. .. GENERATED FROM PYTHON SOURCE LINES 249-258 .. code-block:: Python target_position = SkyCoord(0, 0, unit="deg", frame="galactic") on_radius = 0.2 * u.deg on_region = CircleSkyRegion(center=target_position, radius=on_radius) exclusion_mask = ~geom.to_image().region_mask([on_region]) exclusion_mask.plot() plt.show() .. image-sg:: /tutorials/analysis-3d/images/sphx_glr_cta_data_analysis_003.png :alt: cta data analysis :srcset: /tutorials/analysis-3d/images/sphx_glr_cta_data_analysis_003.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 259-260 Configure spectral analysis .. GENERATED FROM PYTHON SOURCE LINES 260-275 .. code-block:: Python energy_axis = MapAxis.from_energy_bounds(0.1, 40, 40, unit="TeV", name="energy") energy_axis_true = MapAxis.from_energy_bounds( 0.05, 100, 200, unit="TeV", name="energy_true" ) geom = RegionGeom.create(region=on_region, axes=[energy_axis]) dataset_empty = SpectrumDataset.create(geom=geom, energy_axis_true=energy_axis_true) dataset_maker = SpectrumDatasetMaker( containment_correction=False, selection=["counts", "exposure", "edisp"] ) bkg_maker = ReflectedRegionsBackgroundMaker(exclusion_mask=exclusion_mask) safe_mask_masker = SafeMaskMaker(methods=["aeff-max"], aeff_percent=10) .. GENERATED FROM PYTHON SOURCE LINES 276-277 Run data reduction .. GENERATED FROM PYTHON SOURCE LINES 279-289 .. code-block:: Python datasets = Datasets() for observation in observations: dataset = dataset_maker.run( dataset_empty.copy(name=f"obs-{observation.obs_id}"), observation ) dataset_on_off = bkg_maker.run(dataset, observation) dataset_on_off = safe_mask_masker.run(dataset_on_off, observation) datasets.append(dataset_on_off) .. GENERATED FROM PYTHON SOURCE LINES 290-291 Plot results .. GENERATED FROM PYTHON SOURCE LINES 291-300 .. code-block:: Python plt.figure(figsize=(8, 6)) ax = dataset_image.counts.smooth("0.03 deg").plot(vmax=8) on_region.to_pixel(ax.wcs).plot(ax=ax, edgecolor="white") plot_spectrum_datasets_off_regions(datasets, ax=ax) plt.show() .. image-sg:: /tutorials/analysis-3d/images/sphx_glr_cta_data_analysis_004.png :alt: cta data analysis :srcset: /tutorials/analysis-3d/images/sphx_glr_cta_data_analysis_004.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none /home/runner/work/gammapy-docs/gammapy-docs/gammapy/.tox/build_docs/lib/python3.9/site-packages/regions/shapes/circle.py:161: UserWarning: Setting the 'color' property will override the edgecolor or facecolor properties. return Circle(xy=xy, radius=radius, **mpl_kwargs) /home/runner/work/gammapy-docs/gammapy-docs/gammapy/.tox/build_docs/lib/python3.9/site-packages/gammapy/visualization/datasets.py:84: UserWarning: Setting the 'color' property will override the edgecolor or facecolor properties. handle = Patch(**plot_kwargs) .. GENERATED FROM PYTHON SOURCE LINES 301-307 Model fit ~~~~~~~~~ The next step is to fit a spectral model, using all data (i.e. a “global” fit, using all energies). .. GENERATED FROM PYTHON SOURCE LINES 309-322 .. code-block:: Python spectral_model = PowerLawSpectralModel( index=2, amplitude=1e-11 * u.Unit("cm-2 s-1 TeV-1"), reference=1 * u.TeV ) model = SkyModel(spectral_model=spectral_model, name="source-gc") datasets.models = model fit = Fit() result = fit.run(datasets=datasets) print(result) .. rst-class:: sphx-glr-script-out .. code-block:: none OptimizeResult backend : minuit method : migrad success : True message : Optimization terminated successfully. nfev : 104 total stat : 88.36 CovarianceResult backend : minuit method : hesse success : True message : Hesse terminated successfully. .. GENERATED FROM PYTHON SOURCE LINES 323-329 Here we can plot the predicted number of counts for each model and for the background in the dataset. This is especially useful when studying complex field with a lot a sources. There is a function in the visualization sub-package of gammapy that does this automatically. First we need to stack our datasets. .. GENERATED FROM PYTHON SOURCE LINES 329-337 .. code-block:: Python stacked_dataset = datasets.stack_reduce(name="stacked") stacked_dataset.models = model print(stacked_dataset) .. rst-class:: sphx-glr-script-out .. code-block:: none SpectrumDatasetOnOff -------------------- Name : stacked Total counts : 413 Total background counts : 85.43 Total excess counts : 327.57 Predicted counts : 413.95 Predicted background counts : 85.42 Predicted excess counts : 328.53 Exposure min : 9.94e+07 m2 s Exposure max : 2.46e+10 m2 s Number of total bins : 40 Number of fit bins : 30 Fit statistic type : wstat Fit statistic value (-2 log(L)) : 34.70 Number of models : 1 Number of parameters : 3 Number of free parameters : 2 Component 0: SkyModel Name : source-gc Datasets names : None Spectral model type : PowerLawSpectralModel Spatial model type : Temporal model type : Parameters: index : 2.403 +/- 0.06 amplitude : 3.28e-12 +/- 2.3e-13 1 / (cm2 s TeV) reference (frozen): 1.000 TeV Total counts_off : 2095 Acceptance : 88 Acceptance off : 2197 .. GENERATED FROM PYTHON SOURCE LINES 338-340 Call `~gammapy.visualization.plot_npred_signal` to plot the predicted counts. .. GENERATED FROM PYTHON SOURCE LINES 340-346 .. code-block:: Python plot_npred_signal(stacked_dataset) plt.show() .. image-sg:: /tutorials/analysis-3d/images/sphx_glr_cta_data_analysis_005.png :alt: cta data analysis :srcset: /tutorials/analysis-3d/images/sphx_glr_cta_data_analysis_005.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 347-354 Spectral points ~~~~~~~~~~~~~~~ Finally, let’s compute spectral points. The method used is to first choose an energy binning, and then to do a 1-dim likelihood fit / profile to compute the flux and flux error. .. GENERATED FROM PYTHON SOURCE LINES 354-364 .. code-block:: Python # Flux points are computed on stacked datasets energy_edges = MapAxis.from_energy_bounds("1 TeV", "30 TeV", nbin=5).edges fpe = FluxPointsEstimator(energy_edges=energy_edges, source="source-gc") flux_points = fpe.run(datasets=[stacked_dataset]) flux_points.to_table(sed_type="dnde", formatted=True) .. raw:: html
Table length=5
e_refe_mine_maxdndednde_errtssqrt_tsnprednpred_excessstatstat_nullcountssuccess
TeVTeVTeV1 / (cm2 s TeV)1 / (cm2 s TeV)
float64float64float64float64float64float64float64float64[1]float32[1]float64float64float64[1]bool
1.3750.9462.0001.447e-121.783e-13152.51312.350105.7752246942468283.8989213.412165.925106.0True
2.6992.0003.6413.563e-134.835e-14150.65412.27473.0251191207999763.132472.245152.89873.0True
5.2953.6417.7007.332e-141.138e-14121.57011.02653.98359211196148547.455870.624122.19354.0True
11.1987.70016.2846.353e-152.154e-1521.7894.66813.18842984749525210.6604475.74427.53213.0True
21.97116.28429.6451.109e-156.938e-166.2502.5004.14531053887243.1979892.8999.1494.0True


.. GENERATED FROM PYTHON SOURCE LINES 365-372 Plot ~~~~ Let’s plot the spectral model and points. You could do it directly, but for convenience we bundle the model and the flux points in a `~gammapy.datasets.FluxPointsDataset`: .. GENERATED FROM PYTHON SOURCE LINES 372-378 .. code-block:: Python flux_points_dataset = FluxPointsDataset(data=flux_points, models=model) flux_points_dataset.plot_fit() plt.show() .. image-sg:: /tutorials/analysis-3d/images/sphx_glr_cta_data_analysis_006.png :alt: cta data analysis :srcset: /tutorials/analysis-3d/images/sphx_glr_cta_data_analysis_006.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 379-395 Exercises --------- - Re-run the analysis above, varying some analysis parameters, e.g. - Select a few other observations - Change the energy band for the map - Change the spectral model for the fit - Change the energy binning for the spectral points - Change the target. Make a sky image and spectrum for your favourite source. - If you don’t know any, the Crab nebula is the “hello world!” analysis of gamma-ray astronomy. .. GENERATED FROM PYTHON SOURCE LINES 395-400 .. code-block:: Python # print('hello world') # SkyCoord.from_name('crab') .. GENERATED FROM PYTHON SOURCE LINES 401-408 What next? ---------- - This notebook showed an example of a first CTAO analysis with Gammapy, using simulated 1DC data. - Let us know if you have any questions or issues! .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 17.676 seconds) .. _sphx_glr_download_tutorials_analysis-3d_cta_data_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/v1.3?urlpath=lab/tree/notebooks/1.3/tutorials/analysis-3d/cta_data_analysis.ipynb :alt: Launch binder :width: 150 px .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: cta_data_analysis.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: cta_data_analysis.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: cta_data_analysis.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_