.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "tutorials/analysis-3d/simulate_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_simulate_3d.py: 3D map simulation ================= Simulate a 3D observation of a source with the CTA 1DC response and fit it with the assumed source model. Prerequisites ------------- - Knowledge of 3D extraction and datasets used in gammapy, see for example the :doc:`/tutorials/starting/analysis_1` tutorial. Context ------- To simulate a specific observation, it is not always necessary to simulate the full photon list. For many uses cases, simulating directly a reduced binned dataset is enough: the IRFs reduced in the correct geometry are combined with a source model to predict an actual number of counts per bin. The latter is then used to simulate a reduced dataset using Poisson probability distribution. This can be done to check the feasibility of a measurement (performance / sensitivity study), to test whether fitted parameters really provide a good fit to the data etc. Here we will see how to perform a 3D simulation of a CTA observation, assuming both the spectral and spatial morphology of an observed source. **Objective: simulate a 3D observation of a source with CTA using the CTA 1DC response and fit it with the assumed source model.** Proposed approach ----------------- Here we can’t use the regular observation objects that are connected to a `DataStore`. Instead we will create a fake `~gammapy.data.Observation` that contain some pointing information and the CTA 1DC IRFs (that are loaded with `~gammapy.irf.load_irf_dict_from_file`). Then we will create a `~gammapy.datasets.MapDataset` geometry and create it with the `~gammapy.makers.MapDatasetMaker`. Then we will be able to define a model consisting of a `~gammapy.modeling.models.PowerLawSpectralModel` and a `~gammapy.modeling.models.GaussianSpatialModel`. We will assign it to the dataset and fake the count data. .. GENERATED FROM PYTHON SOURCE LINES 52-55 Imports and versions -------------------- .. GENERATED FROM PYTHON SOURCE LINES 55-77 .. code-block:: Python import numpy as np import astropy.units as u from astropy.coordinates import SkyCoord import matplotlib.pyplot as plt # %matplotlib inline from IPython.display import display from gammapy.data import FixedPointingInfo, Observation, observatory_locations from gammapy.datasets import MapDataset from gammapy.irf import load_irf_dict_from_file from gammapy.makers import MapDatasetMaker, SafeMaskMaker from gammapy.maps import MapAxis, WcsGeom from gammapy.modeling import Fit from gammapy.modeling.models import ( FoVBackgroundModel, GaussianSpatialModel, Models, PowerLawSpectralModel, SkyModel, ) .. GENERATED FROM PYTHON SOURCE LINES 78-81 Simulation ---------- .. GENERATED FROM PYTHON SOURCE LINES 84-86 We will simulate using the CTA-1DC IRFs shipped with gammapy .. GENERATED FROM PYTHON SOURCE LINES 86-139 .. code-block:: Python # Loading IRFs irfs = load_irf_dict_from_file( "$GAMMAPY_DATA/cta-1dc/caldb/data/cta/1dc/bcf/South_z20_50h/irf_file.fits" ) # Define the observation parameters (typically the observation duration and the pointing position): livetime = 2.0 * u.hr pointing_position = SkyCoord(0, 0, unit="deg", frame="galactic") # We want to simulate an observation pointing at a fixed position in the sky. # For this, we use the `FixedPointingInfo` class pointing = FixedPointingInfo( fixed_icrs=pointing_position.icrs, ) # Define map geometry for binned simulation energy_reco = MapAxis.from_edges( np.logspace(-1.0, 1.0, 10), unit="TeV", name="energy", interp="log" ) geom = WcsGeom.create( skydir=(0, 0), binsz=0.02, width=(6, 6), frame="galactic", axes=[energy_reco], ) # It is usually useful to have a separate binning for the true energy axis energy_true = MapAxis.from_edges( np.logspace(-1.5, 1.5, 30), unit="TeV", name="energy_true", interp="log" ) empty = MapDataset.create(geom, name="dataset-simu", energy_axis_true=energy_true) # Define sky model to used simulate the data. # Here we use a Gaussian spatial model and a Power Law spectral model. spatial_model = GaussianSpatialModel( lon_0="0.2 deg", lat_0="0.1 deg", sigma="0.3 deg", frame="galactic" ) spectral_model = PowerLawSpectralModel( index=3, amplitude="1e-11 cm-2 s-1 TeV-1", reference="1 TeV" ) model_simu = SkyModel( spatial_model=spatial_model, spectral_model=spectral_model, name="model-simu", ) bkg_model = FoVBackgroundModel(dataset_name="dataset-simu") models = Models([model_simu, bkg_model]) print(models) .. 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) Models Component 0: SkyModel Name : model-simu Datasets names : None Spectral model type : PowerLawSpectralModel Spatial model type : GaussianSpatialModel Temporal model type : Parameters: index : 3.000 +/- 0.00 amplitude : 1.00e-11 +/- 0.0e+00 1 / (cm2 s TeV) reference (frozen): 1.000 TeV lon_0 : 0.200 +/- 0.00 deg lat_0 : 0.100 +/- 0.00 deg sigma : 0.300 +/- 0.00 deg e (frozen): 0.000 phi (frozen): 0.000 deg Component 1: FoVBackgroundModel Name : dataset-simu-bkg Datasets names : ['dataset-simu'] Spectral model type : PowerLawNormSpectralModel Parameters: norm : 1.000 +/- 0.00 tilt (frozen): 0.000 reference (frozen): 1.000 TeV .. GENERATED FROM PYTHON SOURCE LINES 140-146 Now, comes the main part of dataset simulation. We create an in-memory observation and an empty dataset. We then predict the number of counts for the given model, and Poisson fluctuate it using `fake()` to make a simulated counts maps. Keep in mind that it is important to specify the `selection` of the maps that you want to produce .. GENERATED FROM PYTHON SOURCE LINES 146-170 .. code-block:: Python # Create an in-memory observation location = observatory_locations["cta_south"] obs = Observation.create( pointing=pointing, livetime=livetime, irfs=irfs, location=location ) print(obs) # Make the MapDataset maker = MapDatasetMaker(selection=["exposure", "background", "psf", "edisp"]) maker_safe_mask = SafeMaskMaker(methods=["offset-max"], offset_max=4.0 * u.deg) dataset = maker.run(empty, obs) dataset = maker_safe_mask.run(dataset, obs) print(dataset) # Add the model on the dataset and Poisson fluctuate dataset.models = models dataset.fake() # Do a print on the dataset - there is now a counts maps print(dataset) .. rst-class:: sphx-glr-script-out .. code-block:: none Observation obs id : 0 tstart : 51544.00 tstop : 51544.08 duration : 7200.00 s pointing (icrs) : 266.4 deg, -28.9 deg deadtime fraction : 0.0% MapDataset ---------- Name : dataset-simu Total counts : 0 Total background counts : 161250.95 Total excess counts : -161250.95 Predicted counts : 161250.95 Predicted background counts : 161250.95 Predicted excess counts : nan Exposure min : 4.08e+02 m2 s Exposure max : 3.58e+10 m2 s Number of total bins : 810000 Number of fit bins : 804492 Fit statistic type : cash Fit statistic value (-2 log(L)) : nan Number of models : 0 Number of parameters : 0 Number of free parameters : 0 MapDataset ---------- Name : dataset-simu Total counts : 170544 Total background counts : 161250.95 Total excess counts : 9293.05 Predicted counts : 169871.61 Predicted background counts : 161250.95 Predicted excess counts : 8620.66 Exposure min : 4.08e+02 m2 s Exposure max : 3.58e+10 m2 s Number of total bins : 810000 Number of fit bins : 804492 Fit statistic type : cash Fit statistic value (-2 log(L)) : 563788.28 Number of models : 2 Number of parameters : 11 Number of free parameters : 6 Component 0: SkyModel Name : model-simu Datasets names : None Spectral model type : PowerLawSpectralModel Spatial model type : GaussianSpatialModel Temporal model type : Parameters: index : 3.000 +/- 0.00 amplitude : 1.00e-11 +/- 0.0e+00 1 / (cm2 s TeV) reference (frozen): 1.000 TeV lon_0 : 0.200 +/- 0.00 deg lat_0 : 0.100 +/- 0.00 deg sigma : 0.300 +/- 0.00 deg e (frozen): 0.000 phi (frozen): 0.000 deg Component 1: FoVBackgroundModel Name : dataset-simu-bkg Datasets names : ['dataset-simu'] Spectral model type : PowerLawNormSpectralModel Parameters: norm : 1.000 +/- 0.00 tilt (frozen): 0.000 reference (frozen): 1.000 TeV .. GENERATED FROM PYTHON SOURCE LINES 171-176 Now use this dataset as you would in all standard analysis. You can plot the maps, or proceed with your custom analysis. In the next section, we show the standard 3D fitting as in :doc:`/tutorials/analysis-3d/analysis_3d` tutorial. .. GENERATED FROM PYTHON SOURCE LINES 176-182 .. code-block:: Python # To plot, eg, counts: dataset.counts.smooth(0.05 * u.deg).plot_interactive(add_cbar=True, stretch="linear") plt.show() .. image-sg:: /tutorials/analysis-3d/images/sphx_glr_simulate_3d_001.png :alt: simulate 3d :srcset: /tutorials/analysis-3d/images/sphx_glr_simulate_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 - 167 GeV', '167 GeV - 278 GeV', '278 GeV - 464 GeV', '464 GeV - 774 GeV', '774 GeV - 1.29 TeV', '1.29 TeV - 2.15 TeV', '2.15 TeV - 3.59 TeV', '3.59 TeV - 5.99 TeV', '5.99 TeV - 10.0 TeV'), style=SliderStyle(description_width='initial'), value='100 GeV - 167 GeV'), RadioButtons(description='Select stretch:', options=('linear', 'sqrt', 'log'), style=DescriptionStyle(description_width='initial'), value='linear'), Output()), _dom_classes=('widget-interact',)) .. GENERATED FROM PYTHON SOURCE LINES 183-191 Fit --- In this section, we do a usual 3D fit with the same model used to simulated the data and see the stability of the simulations. Often, it is useful to simulate many such datasets and look at the distribution of the reconstructed parameters. .. GENERATED FROM PYTHON SOURCE LINES 191-201 .. code-block:: Python models_fit = models.copy() # We do not want to fit the background in this case, so we will freeze the parameters models_fit["dataset-simu-bkg"].spectral_model.norm.frozen = True models_fit["dataset-simu-bkg"].spectral_model.tilt.frozen = True dataset.models = models_fit print(dataset.models) .. rst-class:: sphx-glr-script-out .. code-block:: none DatasetModels Component 0: SkyModel Name : model-simu Datasets names : None Spectral model type : PowerLawSpectralModel Spatial model type : GaussianSpatialModel Temporal model type : Parameters: index : 3.000 +/- 0.00 amplitude : 1.00e-11 +/- 0.0e+00 1 / (cm2 s TeV) reference (frozen): 1.000 TeV lon_0 : 0.200 +/- 0.00 deg lat_0 : 0.100 +/- 0.00 deg sigma : 0.300 +/- 0.00 deg e (frozen): 0.000 phi (frozen): 0.000 deg Component 1: FoVBackgroundModel Name : dataset-simu-bkg Datasets names : ['dataset-simu'] Spectral model type : PowerLawNormSpectralModel Parameters: norm (frozen): 1.000 tilt (frozen): 0.000 reference (frozen): 1.000 TeV .. GENERATED FROM PYTHON SOURCE LINES 202-209 .. code-block:: Python fit = Fit(optimize_opts={"print_level": 1}) result = fit.run(datasets=[dataset]) dataset.plot_residuals_spatial(method="diff/sqrt(model)", vmin=-0.5, vmax=0.5) plt.show() .. image-sg:: /tutorials/analysis-3d/images/sphx_glr_simulate_3d_002.png :alt: simulate 3d :srcset: /tutorials/analysis-3d/images/sphx_glr_simulate_3d_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 210-212 Compare the injected and fitted models: .. GENERATED FROM PYTHON SOURCE LINES 212-221 .. code-block:: Python print( "True model: \n", model_simu, "\n\n Fitted model: \n", models_fit["model-simu"], ) .. rst-class:: sphx-glr-script-out .. code-block:: none True model: SkyModel Name : model-simu Datasets names : None Spectral model type : PowerLawSpectralModel Spatial model type : GaussianSpatialModel Temporal model type : Parameters: index : 3.000 +/- 0.00 amplitude : 1.00e-11 +/- 0.0e+00 1 / (cm2 s TeV) reference (frozen): 1.000 TeV lon_0 : 0.200 +/- 0.00 deg lat_0 : 0.100 +/- 0.00 deg sigma : 0.300 +/- 0.00 deg e (frozen): 0.000 phi (frozen): 0.000 deg Fitted model: SkyModel Name : model-simu Datasets names : None Spectral model type : PowerLawSpectralModel Spatial model type : GaussianSpatialModel Temporal model type : Parameters: index : 3.007 +/- 0.02 amplitude : 1.01e-11 +/- 3.3e-13 1 / (cm2 s TeV) reference (frozen): 1.000 TeV lon_0 : 0.195 +/- 0.01 deg lat_0 : 0.107 +/- 0.01 deg sigma : 0.301 +/- 0.00 deg e (frozen): 0.000 phi (frozen): 0.000 deg .. GENERATED FROM PYTHON SOURCE LINES 222-224 Get the errors on the fitted parameters from the parameter table .. GENERATED FROM PYTHON SOURCE LINES 224-226 .. code-block:: Python display(result.parameters.to_table()) .. rst-class:: sphx-glr-script-out .. code-block:: none type name value unit error ... frozen is_norm link prior ---- --------- ---------- -------------- --------- ... ------ ------- ---- ----- index 3.0075e+00 1.942e-02 ... False False amplitude 1.0109e-11 cm-2 s-1 TeV-1 3.317e-13 ... False True reference 1.0000e+00 TeV 0.000e+00 ... True False lon_0 1.9477e-01 deg 5.793e-03 ... False False lat_0 1.0717e-01 deg 5.777e-03 ... False False sigma 3.0065e-01 deg 3.948e-03 ... False False e 0.0000e+00 0.000e+00 ... True False phi 0.0000e+00 deg 0.000e+00 ... True False norm 1.0000e+00 0.000e+00 ... True True tilt 0.0000e+00 0.000e+00 ... True False reference 1.0000e+00 TeV 0.000e+00 ... True False .. _sphx_glr_download_tutorials_analysis-3d_simulate_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/simulate_3d.ipynb :alt: Launch binder :width: 150 px .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: simulate_3d.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: simulate_3d.py ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_