.. 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 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, ) from gammapy.visualization import plot_distribution .. GENERATED FROM PYTHON SOURCE LINES 39-46 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 46-84 .. 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/2.0/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 metadata: creator: Gammapy 2.0 date: '2025-08-26T09:05:01.483254' origin: null 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/2.0/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 metadata: creator: Gammapy 2.0 date: '2025-08-26T09:05:01.489044' origin: null .. GENERATED FROM PYTHON SOURCE LINES 85-93 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 93-107 .. code-block:: Python config_stack = config.model_copy(deep=True) config_stack.datasets.stack = True config_joint = config.model_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 108-117 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 117-130 .. code-block:: Python # Reading yaml file: config_stacked = AnalysisConfig.read(path=path / "config_stack.yaml") analysis_stacked = Analysis(config_stacked) # select observations: 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.11/site-packages/astropy/units/core.py:2085: 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.11/site-packages/astropy/units/core.py:2085: 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.11/site-packages/astropy/units/core.py:2085: 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 131-133 We have one final dataset, which we can print and explore .. GENERATED FROM PYTHON SOURCE LINES 133-137 .. 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 138-140 To plot a smooth counts map .. GENERATED FROM PYTHON SOURCE LINES 140-144 .. 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 145-147 And the background map .. GENERATED FROM PYTHON SOURCE LINES 147-151 .. 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 152-154 We can quickly check the PSF .. GENERATED FROM PYTHON SOURCE LINES 154-158 .. 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 159-161 And the energy dispersion in the center of the map .. GENERATED FROM PYTHON SOURCE LINES 161-165 .. 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 166-168 You can also get an excess image with a few lines of code: .. GENERATED FROM PYTHON SOURCE LINES 168-173 .. 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 174-183 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 186-196 To perform the fit on a restricted energy range, we can create a specific *mask*. On the dataset, the `mask_fit` is a `~gammapy.maps.Map` sharing the same geometry as the `~gammapy.datasets.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/details/mask_maps` tutorial. .. GENERATED FROM PYTHON SOURCE LINES 196-228 .. 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 fit = Fit(optimize_opts={"print_level": 1}) result = fit.run(datasets=[dataset_stacked]) .. GENERATED FROM PYTHON SOURCE LINES 229-234 .. _mapdataset_fit_quality: Fit quality assessment and model residuals for a `~gammapy.datasets.MapDataset` ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ .. GENERATED FROM PYTHON SOURCE LINES 237-239 We can access the results dictionary to see if the fit converged: .. GENERATED FROM PYTHON SOURCE LINES 239-243 .. 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 : 189 total stat : 180458.61 CovarianceResult backend : minuit method : hesse success : True message : Hesse terminated successfully. .. GENERATED FROM PYTHON SOURCE LINES 244-246 Check best-fit parameters and error estimates: .. GENERATED FROM PYTHON SOURCE LINES 246-250 .. code-block:: Python display(models_stacked.to_parameters_table()) .. rst-class:: sphx-glr-script-out .. code-block:: none model type name value ... max frozen link prior ----------- ---- --------- ----------- ... --------- ------ ---- ----- gc-source index 2.4160e+00 ... nan False gc-source amplitude 2.6564e-12 ... nan False gc-source reference 1.0000e+00 ... nan True gc-source lambda_ -1.4096e-02 ... nan False gc-source alpha 1.0000e+00 ... nan True gc-source lon_0 -4.8085e-02 ... nan False gc-source lat_0 -5.2598e-02 ... 9.000e+01 False stacked-bkg tilt 0.0000e+00 ... nan True stacked-bkg norm 1.3481e+00 ... nan False stacked-bkg reference 1.0000e+00 ... nan True .. GENERATED FROM PYTHON SOURCE LINES 251-257 A quick way to inspect the model residuals is using the function `~gammapy.datasets.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 257-261 .. 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 262-265 The more general function `~gammapy.datasets.MapDataset.plot_residuals()` can also extract and display spectral residuals in a region: .. GENERATED FROM PYTHON SOURCE LINES 265-274 .. 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 275-286 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 `~gammapy.datasets.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 286-300 .. 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 301-303 Distribution of residuals significance in the full map geometry: .. GENERATED FROM PYTHON SOURCE LINES 303-317 .. 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 318-322 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 325-337 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 337-353 .. 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.11/site-packages/astropy/units/core.py:2085: 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.11/site-packages/astropy/units/core.py:2085: 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.11/site-packages/astropy/units/core.py:2085: 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 : sb1IxzI1 Instrument : CTA Models : Dataset 1: Type : MapDataset Name : N3rVuLt_ Instrument : CTA Models : Dataset 2: Type : MapDataset Name : L77WcQ-7 Instrument : CTA Models : .. GENERATED FROM PYTHON SOURCE LINES 354-356 You can access each one by name or by index, eg: .. GENERATED FROM PYTHON SOURCE LINES 356-360 .. code-block:: Python print(analysis_joint.datasets[0]) .. rst-class:: sphx-glr-script-out .. code-block:: none MapDataset ---------- Name : sb1IxzI1 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 361-366 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 366-387 .. 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 fit_joint = Fit() result_joint = fit_joint.run(datasets=analysis_joint.datasets) .. rst-class:: sphx-glr-script-out .. code-block:: none name counts excess ... n_fit_bins stat_type stat_sum ... -------- ------ ------------------ ... ---------- --------- -------- sb1IxzI1 40481 4466.4930435940405 ... 693940 cash nan N3rVuLt_ 40525 4510.505523195905 ... 693940 cash nan L77WcQ-7 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.416 +/- 0.15 amplitude : 2.66e-12 +/- 3.1e-13 1 / (TeV s cm2) reference (frozen): 1.000 TeV lambda_ : -0.014 +/- 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 : sb1IxzI1-bkg Datasets names : ['sb1IxzI1'] Spectral model type : PowerLawNormSpectralModel Parameters: tilt (frozen): 0.000 norm : 1.000 +/- 0.00 reference (frozen): 1.000 TeV Component 2: FoVBackgroundModel Name : N3rVuLt_-bkg Datasets names : ['N3rVuLt_'] Spectral model type : PowerLawNormSpectralModel Parameters: tilt (frozen): 0.000 norm : 1.000 +/- 0.00 reference (frozen): 1.000 TeV Component 3: FoVBackgroundModel Name : L77WcQ-7-bkg Datasets names : ['L77WcQ-7'] Spectral model type : PowerLawNormSpectralModel Parameters: tilt (frozen): 0.000 norm : 1.000 +/- 0.00 reference (frozen): 1.000 TeV .. GENERATED FROM PYTHON SOURCE LINES 388-393 .. _dataset_fit_quality: Fit quality assessment and model residuals for a joint `~gammapy.datasets.Datasets` ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ .. GENERATED FROM PYTHON SOURCE LINES 396-398 We can access the results dictionary to see if the fit converged: .. GENERATED FROM PYTHON SOURCE LINES 398-402 .. 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 : 225 total stat : 748259.82 CovarianceResult backend : minuit method : hesse success : True message : Hesse terminated successfully. .. GENERATED FROM PYTHON SOURCE LINES 403-405 Check best-fit parameters and error estimates: .. GENERATED FROM PYTHON SOURCE LINES 405-409 .. 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.271 +/- 0.08 amplitude : 2.84e-12 +/- 3.1e-13 1 / (TeV s cm2) 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 : sb1IxzI1-bkg Datasets names : ['sb1IxzI1'] Spectral model type : PowerLawNormSpectralModel Parameters: tilt (frozen): 0.000 norm : 1.118 +/- 0.01 reference (frozen): 1.000 TeV Component 2: FoVBackgroundModel Name : N3rVuLt_-bkg Datasets names : ['N3rVuLt_'] Spectral model type : PowerLawNormSpectralModel Parameters: tilt (frozen): 0.000 norm : 1.119 +/- 0.01 reference (frozen): 1.000 TeV Component 3: FoVBackgroundModel Name : L77WcQ-7-bkg Datasets names : ['L77WcQ-7'] Spectral model type : PowerLawNormSpectralModel Parameters: tilt (frozen): 0.000 norm : 1.111 +/- 0.01 reference (frozen): 1.000 TeV .. GENERATED FROM PYTHON SOURCE LINES 410-416 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 :ref:`mapdataset_fit_quality` - Or, look at a stacked residual map. .. GENERATED FROM PYTHON SOURCE LINES 416-424 .. code-block:: Python stacked = analysis_joint.datasets.stack_reduce() stacked.models = [model_joint] stacked.plot_residuals_spatial(vmin=-1, vmax=1) plt.show() .. 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 .. GENERATED FROM PYTHON SOURCE LINES 425-428 Then, we can access the stacked model residuals as previously shown in the section :ref:`dataset_fit_quality`. .. GENERATED FROM PYTHON SOURCE LINES 431-434 Finally, let us compare the spectral results from the stacked and joint fit: .. GENERATED FROM PYTHON SOURCE LINES 434-452 .. 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 453-462 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 465-468 .. image:: ../../_static/DC1_3d.png .. GENERATED FROM PYTHON SOURCE LINES 471-479 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 28.622 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/v2.0?urlpath=lab/tree/notebooks/2.0/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 ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: analysis_3d.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_