.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "tutorials/analysis-3d/analysis_3d.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_analysis_3d.py: 3D detailed analysis ==================== Perform detailed 3D stacked and joint analysis. This tutorial does a 3D map based analysis on the galactic center, using simulated observations from the CTA-1DC. We will use the high level interface for the data reduction, and then do a detailed modelling. This will be done in two different ways: - stacking all the maps together and fitting the stacked maps - handling all the observations separately and doing a joint fitting on all the maps .. GENERATED FROM PYTHON SOURCE LINES 17-38 .. code-block:: Python from pathlib import Path import numpy as np from scipy.stats import norm import astropy.units as u from regions import CircleSkyRegion # %matplotlib inline import matplotlib.pyplot as plt from IPython.display import display from gammapy.analysis import Analysis, AnalysisConfig from gammapy.estimators import ExcessMapEstimator from gammapy.modeling import Fit from gammapy.modeling.models import ( ExpCutoffPowerLawSpectralModel, FoVBackgroundModel, Models, PointSpatialModel, SkyModel, ) .. GENERATED FROM PYTHON SOURCE LINES 39-41 Check setup ----------- .. GENERATED FROM PYTHON SOURCE LINES 41-47 .. code-block:: Python from gammapy.utils.check import check_tutorials_setup from gammapy.visualization import plot_distribution 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.19 machine : x86_64 system : Linux Gammapy package: version : 1.3.dev414+ge9b22cdff 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.0 astropy : 5.2.2 regions : 0.8 click : 8.1.7 yaml : 6.0.1 IPython : 8.18.1 jupyterlab : not installed matplotlib : 3.9.0 pandas : not installed healpy : 1.16.6 iminuit : 2.25.2 sherpa : 4.16.0 naima : 0.10.0 emcee : 3.1.6 corner : 2.2.2 ray : 2.22.0 Gammapy environment variables: GAMMAPY_DATA : /home/runner/work/gammapy-docs/gammapy-docs/gammapy-datasets/dev .. GENERATED FROM PYTHON SOURCE LINES 48-55 Analysis configuration ---------------------- In this section we select observations and define the analysis geometries, irrespective of joint/stacked analysis. For configuration of the analysis, we will programmatically build a config file from scratch. .. GENERATED FROM PYTHON SOURCE LINES 55-93 .. code-block:: Python config = AnalysisConfig() # The config file is now empty, with only a few defaults specified. print(config) # Selecting the observations config.observations.datastore = "$GAMMAPY_DATA/cta-1dc/index/gps/" config.observations.obs_ids = [110380, 111140, 111159] # Defining a reference geometry for the reduced datasets config.datasets.type = "3d" # Analysis type is 3D config.datasets.geom.wcs.skydir = { "lon": "0 deg", "lat": "0 deg", "frame": "galactic", } # The WCS geometry - centered on the galactic center config.datasets.geom.wcs.width = {"width": "10 deg", "height": "8 deg"} 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.safe_mask.methods = ["aeff-default", "offset-max"] # We now fix the energy axis for the counts map - (the reconstructed energy binning) config.datasets.geom.axes.energy.min = "0.1 TeV" config.datasets.geom.axes.energy.max = "10 TeV" config.datasets.geom.axes.energy.nbins = 10 # We now fix the energy axis for the IRF maps (exposure, etc) - (the true energy binning) config.datasets.geom.axes.energy_true.min = "0.08 TeV" config.datasets.geom.axes.energy_true.max = "12 TeV" config.datasets.geom.axes.energy_true.nbins = 14 print(config) .. rst-class:: sphx-glr-script-out .. code-block:: none AnalysisConfig general: log: {level: info, filename: null, filemode: null, format: null, datefmt: null} outdir: . n_jobs: 1 datasets_file: null models_file: null observations: datastore: /home/runner/work/gammapy-docs/gammapy-docs/gammapy-datasets/dev/hess-dl3-dr1 obs_ids: [] obs_file: null obs_cone: {frame: null, lon: null, lat: null, radius: null} obs_time: {start: null, stop: null} required_irf: [aeff, edisp, psf, bkg] datasets: type: 1d stack: true geom: wcs: skydir: {frame: null, lon: null, lat: null} binsize: 0.02 deg width: {width: 5.0 deg, height: 5.0 deg} binsize_irf: 0.2 deg selection: {offset_max: 2.5 deg} axes: energy: {min: 1.0 TeV, max: 10.0 TeV, nbins: 5} energy_true: {min: 0.5 TeV, max: 20.0 TeV, nbins: 16} map_selection: [counts, exposure, background, psf, edisp] background: method: null exclusion: null parameters: {} safe_mask: methods: [aeff-default] parameters: {} on_region: {frame: null, lon: null, lat: null, radius: null} containment_correction: true fit: fit_range: {min: null, max: null} flux_points: energy: {min: null, max: null, nbins: null} source: source parameters: {selection_optional: all} excess_map: correlation_radius: 0.1 deg parameters: {} energy_edges: {min: null, max: null, nbins: null} light_curve: time_intervals: {start: null, stop: null} energy_edges: {min: null, max: null, nbins: null} source: source parameters: {selection_optional: all} AnalysisConfig general: log: {level: info, filename: null, filemode: null, format: null, datefmt: null} outdir: . n_jobs: 1 datasets_file: null models_file: null observations: datastore: /home/runner/work/gammapy-docs/gammapy-docs/gammapy-datasets/dev/cta-1dc/index/gps obs_ids: [110380, 111140, 111159] obs_file: null obs_cone: {frame: null, lon: null, lat: null, radius: null} obs_time: {start: null, stop: null} required_irf: [aeff, edisp, psf, bkg] datasets: type: 3d stack: true geom: wcs: skydir: {frame: galactic, lon: 0.0 deg, lat: 0.0 deg} binsize: 0.02 deg width: {width: 10.0 deg, height: 8.0 deg} binsize_irf: 0.2 deg selection: {offset_max: 3.5 deg} axes: energy: {min: 0.1 TeV, max: 10.0 TeV, nbins: 10} energy_true: {min: 0.08 TeV, max: 12.0 TeV, nbins: 14} map_selection: [counts, exposure, background, psf, edisp] background: method: null exclusion: null parameters: {} safe_mask: methods: [aeff-default, offset-max] parameters: {} on_region: {frame: null, lon: null, lat: null, radius: null} containment_correction: true fit: fit_range: {min: null, max: null} flux_points: energy: {min: null, max: null, nbins: null} source: source parameters: {selection_optional: all} excess_map: correlation_radius: 0.1 deg parameters: {} energy_edges: {min: null, max: null, nbins: null} light_curve: time_intervals: {start: null, stop: null} energy_edges: {min: null, max: null, nbins: null} source: source parameters: {selection_optional: all} .. GENERATED FROM PYTHON SOURCE LINES 94-102 Configuration for stacked and joint analysis -------------------------------------------- This is done just by specifying the flag on `config.datasets.stack`. Since the internal machinery will work differently for the two cases, we will write it as two config files and save it to disc in YAML format for future reference. .. GENERATED FROM PYTHON SOURCE LINES 102-116 .. code-block:: Python config_stack = config.copy(deep=True) config_stack.datasets.stack = True config_joint = config.copy(deep=True) config_joint.datasets.stack = False # To prevent unnecessary cluttering, we write it in a separate folder. path = Path("analysis_3d") path.mkdir(exist_ok=True) config_joint.write(path=path / "config_joint.yaml", overwrite=True) config_stack.write(path=path / "config_stack.yaml", overwrite=True) .. GENERATED FROM PYTHON SOURCE LINES 117-126 Stacked analysis ---------------- Data reduction ~~~~~~~~~~~~~~ We first show the steps for the stacked analysis and then repeat the same for the joint analysis later .. GENERATED FROM PYTHON SOURCE LINES 126-132 .. code-block:: Python # Reading yaml file: config_stacked = AnalysisConfig.read(path=path / "config_stack.yaml") analysis_stacked = Analysis(config_stacked) .. GENERATED FROM PYTHON SOURCE LINES 133-134 select observations: .. GENERATED FROM PYTHON SOURCE LINES 134-140 .. code-block:: Python analysis_stacked.get_observations() # run data reduction analysis_stacked.get_datasets() .. 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 141-143 We have one final dataset, which we can print and explore .. GENERATED FROM PYTHON SOURCE LINES 143-147 .. code-block:: Python dataset_stacked = analysis_stacked.datasets["stacked"] print(dataset_stacked) .. rst-class:: sphx-glr-script-out .. code-block:: none MapDataset ---------- Name : stacked Total counts : 121241 Total background counts : 108043.52 Total excess counts : 13197.48 Predicted counts : 108043.52 Predicted background counts : 108043.52 Predicted excess counts : nan Exposure min : 6.28e+07 m2 s Exposure max : 1.90e+10 m2 s Number of total bins : 2000000 Number of fit bins : 1411180 Fit statistic type : cash Fit statistic value (-2 log(L)) : nan Number of models : 0 Number of parameters : 0 Number of free parameters : 0 .. GENERATED FROM PYTHON SOURCE LINES 148-150 To plot a smooth counts map .. GENERATED FROM PYTHON SOURCE LINES 150-154 .. code-block:: Python dataset_stacked.counts.smooth(0.02 * u.deg).plot_interactive(add_cbar=True) plt.show() .. image-sg:: /tutorials/analysis-3d/images/sphx_glr_analysis_3d_001.png :alt: analysis 3d :srcset: /tutorials/analysis-3d/images/sphx_glr_analysis_3d_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none interactive(children=(SelectionSlider(continuous_update=False, description='Select energy:', layout=Layout(width='50%'), options=('100 GeV - 158 GeV', '158 GeV - 251 GeV', '251 GeV - 398 GeV', '398 GeV - 631 GeV', '631 GeV - 1.00 TeV', '1.00 TeV - 1.58 TeV', '1.58 TeV - 2.51 TeV', '2.51 TeV - 3.98 TeV', '3.98 TeV - 6.31 TeV', '6.31 TeV - 10.0 TeV'), style=SliderStyle(description_width='initial'), value='100 GeV - 158 GeV'), RadioButtons(description='Select stretch:', index=1, options=('linear', 'sqrt', 'log'), style=DescriptionStyle(description_width='initial'), value='sqrt'), Output()), _dom_classes=('widget-interact',)) .. GENERATED FROM PYTHON SOURCE LINES 155-157 And the background map .. GENERATED FROM PYTHON SOURCE LINES 157-161 .. code-block:: Python dataset_stacked.background.plot_interactive(add_cbar=True) plt.show() .. image-sg:: /tutorials/analysis-3d/images/sphx_glr_analysis_3d_002.png :alt: analysis 3d :srcset: /tutorials/analysis-3d/images/sphx_glr_analysis_3d_002.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none interactive(children=(SelectionSlider(continuous_update=False, description='Select energy:', layout=Layout(width='50%'), options=('100 GeV - 158 GeV', '158 GeV - 251 GeV', '251 GeV - 398 GeV', '398 GeV - 631 GeV', '631 GeV - 1.00 TeV', '1.00 TeV - 1.58 TeV', '1.58 TeV - 2.51 TeV', '2.51 TeV - 3.98 TeV', '3.98 TeV - 6.31 TeV', '6.31 TeV - 10.0 TeV'), style=SliderStyle(description_width='initial'), value='100 GeV - 158 GeV'), RadioButtons(description='Select stretch:', index=1, options=('linear', 'sqrt', 'log'), style=DescriptionStyle(description_width='initial'), value='sqrt'), Output()), _dom_classes=('widget-interact',)) .. GENERATED FROM PYTHON SOURCE LINES 162-164 We can quickly check the PSF .. GENERATED FROM PYTHON SOURCE LINES 164-168 .. code-block:: Python dataset_stacked.psf.peek() plt.show() .. image-sg:: /tutorials/analysis-3d/images/sphx_glr_analysis_3d_003.png :alt: Exposure, Containment radius at 1 TeV, Containment radius at center of map, PSF at center of map :srcset: /tutorials/analysis-3d/images/sphx_glr_analysis_3d_003.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 169-171 And the energy dispersion in the center of the map .. GENERATED FROM PYTHON SOURCE LINES 171-175 .. code-block:: Python dataset_stacked.edisp.peek() plt.show() .. image-sg:: /tutorials/analysis-3d/images/sphx_glr_analysis_3d_004.png :alt: analysis 3d :srcset: /tutorials/analysis-3d/images/sphx_glr_analysis_3d_004.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 176-178 You can also get an excess image with a few lines of code: .. GENERATED FROM PYTHON SOURCE LINES 178-183 .. code-block:: Python excess = dataset_stacked.excess.sum_over_axes() excess.smooth("0.06 deg").plot(stretch="sqrt", add_cbar=True) plt.show() .. image-sg:: /tutorials/analysis-3d/images/sphx_glr_analysis_3d_005.png :alt: analysis 3d :srcset: /tutorials/analysis-3d/images/sphx_glr_analysis_3d_005.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 184-193 Modeling and fitting ~~~~~~~~~~~~~~~~~~~~ Now comes the interesting part of the analysis - choosing appropriate models for our source and fitting them. We choose a point source model with an exponential cutoff power-law spectrum. .. GENERATED FROM PYTHON SOURCE LINES 196-206 To perform the fit on a restricted energy range, we can create a specific *mask*. On the dataset, the `mask_fit` is a `Map` sharing the same geometry as the `MapDataset` and containing boolean data. To create a mask to limit the fit within a restricted energy range, one can rely on the `~gammapy.maps.Geom.energy_mask()` method. For more details on masks and the techniques to create them in gammapy, please checkout the dedicated :doc:`/tutorials/api/mask_maps` tutorial. .. GENERATED FROM PYTHON SOURCE LINES 206-234 .. code-block:: Python dataset_stacked.mask_fit = dataset_stacked.counts.geom.energy_mask( energy_min=0.3 * u.TeV, energy_max=None ) spatial_model = PointSpatialModel( lon_0="-0.05 deg", lat_0="-0.05 deg", frame="galactic" ) spectral_model = ExpCutoffPowerLawSpectralModel( index=2.3, amplitude=2.8e-12 * u.Unit("cm-2 s-1 TeV-1"), reference=1.0 * u.TeV, lambda_=0.02 / u.TeV, ) model = SkyModel( spatial_model=spatial_model, spectral_model=spectral_model, name="gc-source", ) bkg_model = FoVBackgroundModel(dataset_name="stacked") bkg_model.spectral_model.norm.value = 1.3 models_stacked = Models([model, bkg_model]) dataset_stacked.models = models_stacked .. GENERATED FROM PYTHON SOURCE LINES 235-239 .. code-block:: Python fit = Fit(optimize_opts={"print_level": 1}) result = fit.run(datasets=[dataset_stacked]) .. GENERATED FROM PYTHON SOURCE LINES 240-243 Fit quality assessment and model residuals for a `MapDataset` ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ .. GENERATED FROM PYTHON SOURCE LINES 246-248 We can access the results dictionary to see if the fit converged: .. GENERATED FROM PYTHON SOURCE LINES 248-252 .. code-block:: Python print(result) .. rst-class:: sphx-glr-script-out .. code-block:: none OptimizeResult backend : minuit method : migrad success : True message : Optimization terminated successfully. nfev : 183 total stat : 180458.16 CovarianceResult backend : minuit method : hesse success : True message : Hesse terminated successfully. .. GENERATED FROM PYTHON SOURCE LINES 253-255 Check best-fit parameters and error estimates: .. GENERATED FROM PYTHON SOURCE LINES 255-259 .. code-block:: Python display(models_stacked.to_parameters_table()) .. rst-class:: sphx-glr-script-out .. code-block:: none model type name value ... frozen is_norm link prior ----------- ---- --------- ----------- ... ------ ------- ---- ----- gc-source index 2.4144e+00 ... False False gc-source amplitude 2.6623e-12 ... False True gc-source reference 1.0000e+00 ... True False gc-source lambda_ -1.3396e-02 ... False False gc-source alpha 1.0000e+00 ... True False gc-source lon_0 -4.8069e-02 ... False False gc-source lat_0 -5.2610e-02 ... False False stacked-bkg norm 1.3481e+00 ... False True stacked-bkg tilt 0.0000e+00 ... True False stacked-bkg reference 1.0000e+00 ... True False .. GENERATED FROM PYTHON SOURCE LINES 260-266 A quick way to inspect the model residuals is using the function `~MapDataset.plot_residuals_spatial()`. This function computes and plots a residual image (by default, the smoothing radius is `0.1 deg` and `method=diff`, which corresponds to a simple `data - model` plot): .. GENERATED FROM PYTHON SOURCE LINES 266-270 .. code-block:: Python dataset_stacked.plot_residuals_spatial(method="diff/sqrt(model)", vmin=-1, vmax=1) plt.show() .. image-sg:: /tutorials/analysis-3d/images/sphx_glr_analysis_3d_006.png :alt: analysis 3d :srcset: /tutorials/analysis-3d/images/sphx_glr_analysis_3d_006.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 271-274 The more general function `~MapDataset.plot_residuals()` can also extract and display spectral residuals in a region: .. GENERATED FROM PYTHON SOURCE LINES 274-283 .. code-block:: Python region = CircleSkyRegion(spatial_model.position, radius=0.15 * u.deg) dataset_stacked.plot_residuals( kwargs_spatial=dict(method="diff/sqrt(model)", vmin=-1, vmax=1), kwargs_spectral=dict(region=region), ) plt.show() .. image-sg:: /tutorials/analysis-3d/images/sphx_glr_analysis_3d_007.png :alt: analysis 3d :srcset: /tutorials/analysis-3d/images/sphx_glr_analysis_3d_007.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 284-293 This way of accessing residuals is quick and handy, but comes with limitations. For example: - In case a fitting energy range was defined using a `MapDataset.mask_fit`, it won’t be taken into account. Residuals will be summed up over the whole reconstructed energy range - In order to make a proper statistic treatment, instead of simple residuals a proper residuals significance map should be computed A more accurate way to inspect spatial residuals is the following: .. GENERATED FROM PYTHON SOURCE LINES 293-307 .. code-block:: Python estimator = ExcessMapEstimator( correlation_radius="0.1 deg", selection_optional=[], energy_edges=[0.1, 1, 10] * u.TeV, ) result = estimator.run(dataset_stacked) result["sqrt_ts"].plot_grid( figsize=(12, 4), cmap="coolwarm", add_cbar=True, vmin=-5, vmax=5, ncols=2 ) plt.show() .. image-sg:: /tutorials/analysis-3d/images/sphx_glr_analysis_3d_008.png :alt: Energy 100 GeV - 1.00 TeV, Energy 1.00 TeV - 10.0 TeV :srcset: /tutorials/analysis-3d/images/sphx_glr_analysis_3d_008.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 308-310 Distribution of residuals significance in the full map geometry: .. GENERATED FROM PYTHON SOURCE LINES 310-324 .. code-block:: Python significance_map = result["sqrt_ts"] kwargs_hist = {"density": True, "alpha": 0.9, "color": "red", "bins": 40} ax, res = plot_distribution( significance_map, func="norm", kwargs_hist=kwargs_hist, kwargs_axes={"xlim": (-5, 5)}, ) plt.show() .. image-sg:: /tutorials/analysis-3d/images/sphx_glr_analysis_3d_009.png :alt: analysis 3d :srcset: /tutorials/analysis-3d/images/sphx_glr_analysis_3d_009.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 325-329 Here we could also plot the number of predicted counts for each model and for the background in our dataset by using the `~gammapy.visualization.plot_npred_signal` function. .. GENERATED FROM PYTHON SOURCE LINES 332-344 Joint analysis -------------- In this section, we perform a joint analysis of the same data. Of course, joint fitting is considerably heavier than stacked one, and should always be handled with care. For brevity, we only show the analysis for a point source fitting without re-adding a diffuse component again. Data reduction ~~~~~~~~~~~~~~ .. GENERATED FROM PYTHON SOURCE LINES 346-361 .. code-block:: Python # Read the yaml file from disk config_joint = AnalysisConfig.read(path=path / "config_joint.yaml") analysis_joint = Analysis(config_joint) # select observations: analysis_joint.get_observations() # run data reduction analysis_joint.get_datasets() # You can see there are 3 datasets now print(analysis_joint.datasets) .. 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) Datasets -------- Dataset 0: Type : MapDataset Name : j9ENGRHT Instrument : CTA Models : Dataset 1: Type : MapDataset Name : Jocs6eXK Instrument : CTA Models : Dataset 2: Type : MapDataset Name : nKgXgwU2 Instrument : CTA Models : .. GENERATED FROM PYTHON SOURCE LINES 362-364 You can access each one by name or by index, eg: .. GENERATED FROM PYTHON SOURCE LINES 364-368 .. code-block:: Python print(analysis_joint.datasets[0]) .. rst-class:: sphx-glr-script-out .. code-block:: none MapDataset ---------- Name : j9ENGRHT Total counts : 40481 Total background counts : 36014.51 Total excess counts : 4466.49 Predicted counts : 36014.51 Predicted background counts : 36014.51 Predicted excess counts : nan Exposure min : 6.28e+07 m2 s Exposure max : 6.68e+09 m2 s Number of total bins : 1085000 Number of fit bins : 693940 Fit statistic type : cash Fit statistic value (-2 log(L)) : nan Number of models : 0 Number of parameters : 0 Number of free parameters : 0 .. GENERATED FROM PYTHON SOURCE LINES 369-374 After the data reduction stage, it is nice to get a quick summary info on the datasets. Here, we look at the statistics in the center of Map, by passing an appropriate `region`. To get info on the entire spatial map, omit the region argument. .. GENERATED FROM PYTHON SOURCE LINES 374-391 .. code-block:: Python display(analysis_joint.datasets.info_table()) models_joint = Models() model_joint = model.copy(name="source-joint") models_joint.append(model_joint) for dataset in analysis_joint.datasets: bkg_model = FoVBackgroundModel(dataset_name=dataset.name) models_joint.append(bkg_model) print(models_joint) # and set the new model analysis_joint.datasets.models = models_joint .. rst-class:: sphx-glr-script-out .. code-block:: none name counts excess ... n_fit_bins stat_type stat_sum ... -------- ------ ------------------ ... ---------- --------- -------- j9ENGRHT 40481 4466.4930435940405 ... 693940 cash nan Jocs6eXK 40525 4510.505523195905 ... 693940 cash nan nKgXgwU2 40235 4220.480554966147 ... 693940 cash nan Models Component 0: SkyModel Name : source-joint Datasets names : None Spectral model type : ExpCutoffPowerLawSpectralModel Spatial model type : PointSpatialModel Temporal model type : Parameters: index : 2.414 +/- 0.15 amplitude : 2.66e-12 +/- 3.1e-13 1 / (cm2 s TeV) reference (frozen): 1.000 TeV lambda_ : -0.013 +/- 0.07 1 / TeV alpha (frozen): 1.000 lon_0 : -0.048 +/- 0.00 deg lat_0 : -0.053 +/- 0.00 deg Component 1: FoVBackgroundModel Name : j9ENGRHT-bkg Datasets names : ['j9ENGRHT'] Spectral model type : PowerLawNormSpectralModel Parameters: norm : 1.000 +/- 0.00 tilt (frozen): 0.000 reference (frozen): 1.000 TeV Component 2: FoVBackgroundModel Name : Jocs6eXK-bkg Datasets names : ['Jocs6eXK'] Spectral model type : PowerLawNormSpectralModel Parameters: norm : 1.000 +/- 0.00 tilt (frozen): 0.000 reference (frozen): 1.000 TeV Component 3: FoVBackgroundModel Name : nKgXgwU2-bkg Datasets names : ['nKgXgwU2'] Spectral model type : PowerLawNormSpectralModel Parameters: norm : 1.000 +/- 0.00 tilt (frozen): 0.000 reference (frozen): 1.000 TeV .. GENERATED FROM PYTHON SOURCE LINES 392-396 .. code-block:: Python fit_joint = Fit() result_joint = fit_joint.run(datasets=analysis_joint.datasets) .. GENERATED FROM PYTHON SOURCE LINES 397-400 Fit quality assessment and model residuals for a joint `Datasets` ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ .. GENERATED FROM PYTHON SOURCE LINES 403-405 We can access the results dictionary to see if the fit converged: .. GENERATED FROM PYTHON SOURCE LINES 405-409 .. code-block:: Python print(result_joint) .. rst-class:: sphx-glr-script-out .. code-block:: none OptimizeResult backend : minuit method : migrad success : True message : Optimization terminated successfully. nfev : 220 total stat : 748259.32 CovarianceResult backend : minuit method : hesse success : True message : Hesse terminated successfully. .. GENERATED FROM PYTHON SOURCE LINES 410-412 Check best-fit parameters and error estimates: .. GENERATED FROM PYTHON SOURCE LINES 412-416 .. code-block:: Python print(models_joint) .. rst-class:: sphx-glr-script-out .. code-block:: none Models Component 0: SkyModel Name : source-joint Datasets names : None Spectral model type : ExpCutoffPowerLawSpectralModel Spatial model type : PointSpatialModel Temporal model type : Parameters: index : 2.272 +/- 0.08 amplitude : 2.84e-12 +/- 3.1e-13 1 / (cm2 s TeV) reference (frozen): 1.000 TeV lambda_ : 0.039 +/- 0.05 1 / TeV alpha (frozen): 1.000 lon_0 : -0.049 +/- 0.00 deg lat_0 : -0.053 +/- 0.00 deg Component 1: FoVBackgroundModel Name : j9ENGRHT-bkg Datasets names : ['j9ENGRHT'] Spectral model type : PowerLawNormSpectralModel Parameters: norm : 1.118 +/- 0.01 tilt (frozen): 0.000 reference (frozen): 1.000 TeV Component 2: FoVBackgroundModel Name : Jocs6eXK-bkg Datasets names : ['Jocs6eXK'] Spectral model type : PowerLawNormSpectralModel Parameters: norm : 1.119 +/- 0.01 tilt (frozen): 0.000 reference (frozen): 1.000 TeV Component 3: FoVBackgroundModel Name : nKgXgwU2-bkg Datasets names : ['nKgXgwU2'] Spectral model type : PowerLawNormSpectralModel Parameters: norm : 1.111 +/- 0.01 tilt (frozen): 0.000 reference (frozen): 1.000 TeV .. GENERATED FROM PYTHON SOURCE LINES 417-423 Since the joint dataset is made of multiple datasets, we can either: - Look at the residuals for each dataset separately. In this case, we can directly refer to the section `Fit quality and model residuals for a MapDataset` in this notebook - Or, look at a stacked residual map. .. GENERATED FROM PYTHON SOURCE LINES 423-431 .. code-block:: Python stacked = analysis_joint.datasets.stack_reduce() stacked.models = [model_joint] plt.figure() stacked.plot_residuals_spatial(vmin=-1, vmax=1) .. image-sg:: /tutorials/analysis-3d/images/sphx_glr_analysis_3d_010.png :alt: analysis 3d :srcset: /tutorials/analysis-3d/images/sphx_glr_analysis_3d_010.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 432-436 Then, we can access the stacked model residuals as previously shown in the section `Fit quality and model residuals for a MapDataset` in this notebook. .. GENERATED FROM PYTHON SOURCE LINES 439-442 Finally, let us compare the spectral results from the stacked and joint fit: .. GENERATED FROM PYTHON SOURCE LINES 442-460 .. code-block:: Python def plot_spectrum(model, ax, label, color): spec = model.spectral_model energy_bounds = [0.3, 10] * u.TeV spec.plot( ax=ax, energy_bounds=energy_bounds, energy_power=2, label=label, color=color ) spec.plot_error(ax=ax, energy_bounds=energy_bounds, energy_power=2, color=color) fig, ax = plt.subplots() plot_spectrum(model, ax=ax, label="stacked", color="tab:blue") plot_spectrum(model_joint, ax=ax, label="joint", color="tab:orange") ax.legend() plt.show() .. image-sg:: /tutorials/analysis-3d/images/sphx_glr_analysis_3d_011.png :alt: analysis 3d :srcset: /tutorials/analysis-3d/images/sphx_glr_analysis_3d_011.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 461-470 Summary ------- Note that this notebook aims to show you the procedure of a 3D analysis using just a few observations. Results get much better for a more complete analysis considering the GPS dataset from the CTA First Data Challenge (DC-1) and also the CTA model for the Galactic diffuse emission, as shown in the next image: .. GENERATED FROM PYTHON SOURCE LINES 473-475 .. image:: ../../_static/DC1_3d.png .. GENERATED FROM PYTHON SOURCE LINES 478-486 Exercises --------- - Analyse the second source in the field of view: G0.9+0.1 and add it to the combined model. - Perform modeling in more details - Add diffuse component, get flux points. .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 32.046 seconds) .. _sphx_glr_download_tutorials_analysis-3d_analysis_3d.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-3d/analysis_3d.ipynb :alt: Launch binder :width: 150 px .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: analysis_3d.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: analysis_3d.py ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_