.. 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-36 .. 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, ) .. GENERATED FROM PYTHON SOURCE LINES 37-39 Check setup ----------- .. GENERATED FROM PYTHON SOURCE LINES 39-45 .. 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.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 46-53 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 53-91 .. 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/1.3/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 1.3 date: '2024-11-26T10:11:39.316231' 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/1.3/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 1.3 date: '2024-11-26T10:11:39.324287' origin: null .. GENERATED FROM PYTHON SOURCE LINES 92-100 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 100-114 .. 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) .. rst-class:: sphx-glr-script-out .. code-block:: none /home/runner/work/gammapy-docs/gammapy-docs/gammapy/examples/tutorials/analysis-3d/analysis_3d.py:101: PydanticDeprecatedSince20: The `copy` method is deprecated; use `model_copy` instead. See the docstring of `BaseModel.copy` for details about how to handle `include` and `exclude`. Deprecated in Pydantic V2.0 to be removed in V3.0. See Pydantic V2 Migration Guide at https://errors.pydantic.dev/2.10/migration/ config_stack = config.copy(deep=True) /home/runner/work/gammapy-docs/gammapy-docs/gammapy/examples/tutorials/analysis-3d/analysis_3d.py:104: PydanticDeprecatedSince20: The `copy` method is deprecated; use `model_copy` instead. See the docstring of `BaseModel.copy` for details about how to handle `include` and `exclude`. Deprecated in Pydantic V2.0 to be removed in V3.0. See Pydantic V2 Migration Guide at https://errors.pydantic.dev/2.10/migration/ config_joint = config.copy(deep=True) .. GENERATED FROM PYTHON SOURCE LINES 115-124 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 124-130 .. 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 131-132 select observations: .. GENERATED FROM PYTHON SOURCE LINES 132-138 .. 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 139-141 We have one final dataset, which we can print and explore .. GENERATED FROM PYTHON SOURCE LINES 141-145 .. 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 146-148 To plot a smooth counts map .. GENERATED FROM PYTHON SOURCE LINES 148-152 .. 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 153-155 And the background map .. GENERATED FROM PYTHON SOURCE LINES 155-159 .. 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 160-162 We can quickly check the PSF .. GENERATED FROM PYTHON SOURCE LINES 162-166 .. 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 167-169 And the energy dispersion in the center of the map .. GENERATED FROM PYTHON SOURCE LINES 169-173 .. 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 174-176 You can also get an excess image with a few lines of code: .. GENERATED FROM PYTHON SOURCE LINES 176-181 .. 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 182-191 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 194-204 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 204-232 .. 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 233-237 .. code-block:: Python fit = Fit(optimize_opts={"print_level": 1}) result = fit.run(datasets=[dataset_stacked]) .. GENERATED FROM PYTHON SOURCE LINES 238-241 Fit quality assessment and model residuals for a `MapDataset` ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ .. GENERATED FROM PYTHON SOURCE LINES 244-246 We can access the results dictionary to see if the fit converged: .. GENERATED FROM PYTHON SOURCE LINES 246-250 .. 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 : 184 total stat : 180458.61 CovarianceResult backend : minuit method : hesse success : True message : Hesse terminated successfully. .. GENERATED FROM PYTHON SOURCE LINES 251-253 Check best-fit parameters and error estimates: .. GENERATED FROM PYTHON SOURCE LINES 253-257 .. 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.4161e+00 ... nan False gc-source amplitude 2.6561e-12 ... nan False gc-source reference 1.0000e+00 ... nan True gc-source lambda_ -1.4151e-02 ... nan False gc-source alpha 1.0000e+00 ... nan True gc-source lon_0 -4.8084e-02 ... nan False gc-source lat_0 -5.2595e-02 ... 9.000e+01 False stacked-bkg norm 1.3481e+00 ... nan False stacked-bkg tilt 0.0000e+00 ... nan True stacked-bkg reference 1.0000e+00 ... nan True .. GENERATED FROM PYTHON SOURCE LINES 258-264 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 264-268 .. 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 269-272 The more general function `~MapDataset.plot_residuals()` can also extract and display spectral residuals in a region: .. GENERATED FROM PYTHON SOURCE LINES 272-281 .. 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 282-291 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 291-305 .. 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 306-308 Distribution of residuals significance in the full map geometry: .. GENERATED FROM PYTHON SOURCE LINES 308-322 .. 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 323-327 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 330-342 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 344-359 .. 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 : suqll3vI Instrument : CTA Models : Dataset 1: Type : MapDataset Name : Da-Bk9ou Instrument : CTA Models : Dataset 2: Type : MapDataset Name : G6s80VnL Instrument : CTA Models : .. GENERATED FROM PYTHON SOURCE LINES 360-362 You can access each one by name or by index, eg: .. GENERATED FROM PYTHON SOURCE LINES 362-366 .. code-block:: Python print(analysis_joint.datasets[0]) .. rst-class:: sphx-glr-script-out .. code-block:: none MapDataset ---------- Name : suqll3vI 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 367-372 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 372-389 .. 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 ... -------- ------ ----------------- ... ---------- --------- -------- suqll3vI 40481 4466.493043594026 ... 693940 cash nan Da-Bk9ou 40525 4510.505523195905 ... 693940 cash nan G6s80VnL 40235 4220.480554966154 ... 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 / (cm2 s TeV) 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 : suqll3vI-bkg Datasets names : ['suqll3vI'] Spectral model type : PowerLawNormSpectralModel Parameters: norm : 1.000 +/- 0.00 tilt (frozen): 0.000 reference (frozen): 1.000 TeV Component 2: FoVBackgroundModel Name : Da-Bk9ou-bkg Datasets names : ['Da-Bk9ou'] Spectral model type : PowerLawNormSpectralModel Parameters: norm : 1.000 +/- 0.00 tilt (frozen): 0.000 reference (frozen): 1.000 TeV Component 3: FoVBackgroundModel Name : G6s80VnL-bkg Datasets names : ['G6s80VnL'] Spectral model type : PowerLawNormSpectralModel Parameters: norm : 1.000 +/- 0.00 tilt (frozen): 0.000 reference (frozen): 1.000 TeV .. GENERATED FROM PYTHON SOURCE LINES 390-394 .. code-block:: Python fit_joint = Fit() result_joint = fit_joint.run(datasets=analysis_joint.datasets) .. GENERATED FROM PYTHON SOURCE LINES 395-398 Fit quality assessment and model residuals for a joint `Datasets` ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ .. GENERATED FROM PYTHON SOURCE LINES 401-403 We can access the results dictionary to see if the fit converged: .. GENERATED FROM PYTHON SOURCE LINES 403-407 .. 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 : 238 total stat : 748259.82 CovarianceResult backend : minuit method : hesse success : True message : Hesse terminated successfully. .. GENERATED FROM PYTHON SOURCE LINES 408-410 Check best-fit parameters and error estimates: .. GENERATED FROM PYTHON SOURCE LINES 410-414 .. 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 / (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 : suqll3vI-bkg Datasets names : ['suqll3vI'] Spectral model type : PowerLawNormSpectralModel Parameters: norm : 1.118 +/- 0.01 tilt (frozen): 0.000 reference (frozen): 1.000 TeV Component 2: FoVBackgroundModel Name : Da-Bk9ou-bkg Datasets names : ['Da-Bk9ou'] Spectral model type : PowerLawNormSpectralModel Parameters: norm : 1.119 +/- 0.01 tilt (frozen): 0.000 reference (frozen): 1.000 TeV Component 3: FoVBackgroundModel Name : G6s80VnL-bkg Datasets names : ['G6s80VnL'] Spectral model type : PowerLawNormSpectralModel Parameters: norm : 1.111 +/- 0.01 tilt (frozen): 0.000 reference (frozen): 1.000 TeV .. GENERATED FROM PYTHON SOURCE LINES 415-421 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 421-429 .. 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 430-434 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 437-440 Finally, let us compare the spectral results from the stacked and joint fit: .. GENERATED FROM PYTHON SOURCE LINES 440-458 .. 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 459-468 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 471-473 .. image:: ../../_static/DC1_3d.png .. GENERATED FROM PYTHON SOURCE LINES 476-484 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 35.506 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/v1.3?urlpath=lab/tree/notebooks/1.3/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 `_