.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "tutorials/analysis-time/light_curve_simulation.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-time_light_curve_simulation.py: Simulating and fitting a time varying source ============================================ Simulate and fit a time decaying light curve of a source using the CTA 1DC response. Prerequisites ------------- - To understand how a single binned simulation works, please refer to :doc:`/tutorials/analysis-1d/spectrum_simulation` tutorial and :doc:`/tutorials/analysis-3d/simulate_3d` tutorial for 1D and 3D simulations, respectively. - For details of light curve extraction using gammapy, refer to the two tutorials :doc:`/tutorials/analysis-time/light_curve` and :doc:`/tutorials/analysis-time/light_curve_flare`. Context ------- Frequently, studies of variable sources (eg: decaying GRB light curves, AGN flares, etc.) require time variable simulations. For most use cases, generating an event list is an overkill, and it suffices to use binned simulations using a temporal model. **Objective: Simulate and fit a time decaying light curve of a source with CTA using the CTA 1DC response.** Proposed approach ----------------- We will simulate 10 spectral datasets within given time intervals (Good Time Intervals) following a given spectral (a power law) and temporal profile (an exponential decay, with a decay time of 6 hr). These are then analysed using the light curve estimator to obtain flux points. Modelling and fitting of lightcurves can be done either - directly on the output of the `~gammapy.estimators.LightCurveEstimator` (at the DL5 level) - fit the simulated datasets (at the DL4 level) In summary, the necessary steps are: - Choose observation parameters including a list of `gammapy.data.GTI` - Define temporal and spectral models from the :ref:`model-gallery` as per science case - Perform the simulation (in 1D or 3D) - Extract the light curve from the reduced dataset as shown in :doc:`/tutorials/analysis-time/light_curve` tutorial - Optionally, we show here how to fit the simulated datasets using a source model Setup ----- As usual, we’ll start with some general imports… .. GENERATED FROM PYTHON SOURCE LINES 59-73 .. code-block:: Python import logging import numpy as np import astropy.units as u from astropy.coordinates import SkyCoord from astropy.time import Time # %matplotlib inline import matplotlib.pyplot as plt from IPython.display import display log = logging.getLogger(__name__) .. GENERATED FROM PYTHON SOURCE LINES 74-76 And some gammapy specific imports .. GENERATED FROM PYTHON SOURCE LINES 76-95 .. code-block:: Python import warnings from gammapy.data import FixedPointingInfo, Observation, observatory_locations from gammapy.datasets import Datasets, FluxPointsDataset, SpectrumDataset from gammapy.estimators import LightCurveEstimator from gammapy.irf import load_irf_dict_from_file from gammapy.makers import SpectrumDatasetMaker from gammapy.maps import MapAxis, RegionGeom, TimeMapAxis from gammapy.modeling import Fit from gammapy.modeling.models import ( ExpDecayTemporalModel, PowerLawSpectralModel, SkyModel, ) warnings.filterwarnings( action="ignore", message="overflow encountered in exp", module="astropy" ) .. GENERATED FROM PYTHON SOURCE LINES 96-98 Check setup ----------- .. GENERATED FROM PYTHON SOURCE LINES 98-103 .. code-block:: Python from gammapy.utils.check import check_tutorials_setup 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 104-106 We first define our preferred time format: .. GENERATED FROM PYTHON SOURCE LINES 106-110 .. code-block:: Python TimeMapAxis.time_format = "iso" .. GENERATED FROM PYTHON SOURCE LINES 111-120 Simulating a light curve ------------------------ We will simulate 10 spectra between 300 GeV and 10 TeV using an `~gammapy.modeling.models.PowerLawSpectralModel` and a `~gammapy.modeling.models.ExpDecayTemporalModel`. The important thing to note here is how to attach a different `GTI` to each dataset. Since we use spectrum datasets here, we will use a `~gammapy.maps.RegionGeom`. .. GENERATED FROM PYTHON SOURCE LINES 120-142 .. 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" ) # Reconstructed and true energy axis energy_axis = MapAxis.from_edges( np.logspace(-0.5, 1.0, 10), unit="TeV", name="energy", interp="log" ) energy_axis_true = MapAxis.from_edges( np.logspace(-1.2, 2.0, 31), unit="TeV", name="energy_true", interp="log" ) geom = RegionGeom.create("galactic;circle(0, 0, 0.11)", axes=[energy_axis]) # Pointing position to be supplied as a `FixedPointingInfo` pointing = FixedPointingInfo( fixed_icrs=SkyCoord(0.5, 0.5, unit="deg", frame="galactic").icrs, ) .. 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) .. GENERATED FROM PYTHON SOURCE LINES 143-147 Note that observations are usually conducted in Wobble mode, in which the source is not in the center of the camera. This allows to have a symmetrical sky position from which background can be estimated. .. GENERATED FROM PYTHON SOURCE LINES 147-166 .. code-block:: Python # Define the source model: A combination of spectral and temporal model gti_t0 = Time("2020-03-01") spectral_model = PowerLawSpectralModel( index=3, amplitude="1e-11 cm-2 s-1 TeV-1", reference="1 TeV" ) temporal_model = ExpDecayTemporalModel(t0="6 h", t_ref=gti_t0.mjd * u.d) model_simu = SkyModel( spectral_model=spectral_model, temporal_model=temporal_model, name="model-simu", ) # Look at the model display(model_simu.parameters.to_table()) .. rst-class:: sphx-glr-script-out .. code-block:: none type name value unit error ... frozen is_norm link prior ---- --------- ---------- -------------- --------- ... ------ ------- ---- ----- index 3.0000e+00 0.000e+00 ... False False amplitude 1.0000e-11 cm-2 s-1 TeV-1 0.000e+00 ... False True reference 1.0000e+00 TeV 0.000e+00 ... True False t0 6.0000e+00 h 0.000e+00 ... False False t_ref 5.8909e+04 d 0.000e+00 ... True False .. GENERATED FROM PYTHON SOURCE LINES 167-170 Now, define the start and observation livetime wrt to the reference time, ``gti_t0`` .. GENERATED FROM PYTHON SOURCE LINES 170-177 .. code-block:: Python n_obs = 10 tstart = gti_t0 + [1, 2, 3, 5, 8, 10, 20, 22, 23, 24] * u.h lvtm = [55, 25, 26, 40, 40, 50, 40, 52, 43, 47] * u.min .. GENERATED FROM PYTHON SOURCE LINES 178-180 Now perform the simulations .. GENERATED FROM PYTHON SOURCE LINES 180-207 .. code-block:: Python datasets = Datasets() empty = SpectrumDataset.create( geom=geom, energy_axis_true=energy_axis_true, name="empty" ) maker = SpectrumDatasetMaker(selection=["exposure", "background", "edisp"]) for idx in range(n_obs): obs = Observation.create( pointing=pointing, livetime=lvtm[idx], tstart=tstart[idx], irfs=irfs, reference_time=gti_t0, obs_id=idx, location=observatory_locations["cta_south"], ) empty_i = empty.copy(name=f"dataset-{idx}") dataset = maker.run(empty_i, obs) dataset.models = model_simu dataset.fake() datasets.append(dataset) .. GENERATED FROM PYTHON SOURCE LINES 208-211 The reduced datasets have been successfully simulated. Let’s take a quick look into our datasets. .. GENERATED FROM PYTHON SOURCE LINES 211-215 .. code-block:: Python display(datasets.info_table()) .. rst-class:: sphx-glr-script-out .. code-block:: none name counts excess ... n_fit_bins stat_type stat_sum ... --------- ------ ------------------ ... ---------- --------- ------------------- dataset-0 812 791.6812310184489 ... 9 cash -6544.801116964405 dataset-1 305 295.76419591747674 ... 9 cash -1848.5164045025529 dataset-2 269 259.3947637541758 ... 9 cash -1553.926967307144 dataset-3 301 286.2227134679628 ... 9 cash -1811.1958962266817 dataset-4 195 180.22271346796282 ... 9 cash -1018.6102107826875 dataset-5 183 164.52839183495354 ... 9 cash -912.4778519565482 dataset-6 42 27.22271346796283 ... 9 cash -97.39153688836538 dataset-7 49 29.78952750835169 ... 9 cash -122.55922091763816 dataset-8 24 8.114416978060046 ... 9 cash -30.787758061870708 dataset-9 36 18.636688324856333 ... 9 cash -66.57719432137594 .. GENERATED FROM PYTHON SOURCE LINES 216-225 Extract the lightcurve ---------------------- This section uses standard light curve estimation tools for a 1D extraction. Only a spectral model needs to be defined in this case. Since the estimator returns the integrated flux separately for each time bin, the temporal model need not be accounted for at this stage. We extract the lightcurve in 3 energy bins. .. GENERATED FROM PYTHON SOURCE LINES 225-235 .. code-block:: Python # Define the model: spectral_model = PowerLawSpectralModel( index=3, amplitude="1e-11 cm-2 s-1 TeV-1", reference="1 TeV" ) model_fit = SkyModel(spectral_model=spectral_model, name="model-fit") # Attach model to all datasets datasets.models = model_fit .. GENERATED FROM PYTHON SOURCE LINES 236-251 .. code-block:: Python lc_maker_1d = LightCurveEstimator( energy_edges=[0.3, 0.6, 1.0, 10] * u.TeV, source="model-fit", selection_optional=["ul"], ) lc_1d = lc_maker_1d.run(datasets) fig, ax = plt.subplots( figsize=(8, 6), gridspec_kw={"left": 0.16, "bottom": 0.2, "top": 0.98, "right": 0.98}, ) lc_1d.plot(ax=ax, marker="o", axis_name="time", sed_type="flux") plt.show() .. image-sg:: /tutorials/analysis-time/images/sphx_glr_light_curve_simulation_001.png :alt: light curve simulation :srcset: /tutorials/analysis-time/images/sphx_glr_light_curve_simulation_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 252-263 Fitting temporal models ----------------------- We have the reconstructed lightcurve at this point. Now we want to fit a profile to the obtained light curves, using a joint fitting across the different datasets, while simultaneously minimising across the temporal model parameters as well. The temporal models can be applied - directly on the obtained lightcurve - on the simulated datasets .. GENERATED FROM PYTHON SOURCE LINES 266-271 Fitting the obtained light curve ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ We will first convert the obtained light curve to a `~gammapy.datasets.FluxPointsDataset` and fit it with a spectral and temporal model .. GENERATED FROM PYTHON SOURCE LINES 271-276 .. code-block:: Python # Create the datasets by iterating over the returned lightcurve dataset_fp = FluxPointsDataset(data=lc_1d, name="dataset_lc") .. GENERATED FROM PYTHON SOURCE LINES 277-281 We will fit the amplitude, spectral index and the decay time scale. Note that ``t_ref`` should be fixed by default for the `~gammapy.modeling.models.ExpDecayTemporalModel`. .. GENERATED FROM PYTHON SOURCE LINES 281-299 .. code-block:: Python # Define the model: spectral_model1 = PowerLawSpectralModel( index=2.0, amplitude="1e-12 cm-2 s-1 TeV-1", reference="1 TeV" ) temporal_model1 = ExpDecayTemporalModel(t0="10 h", t_ref=gti_t0.mjd * u.d) model = SkyModel( spectral_model=spectral_model1, temporal_model=temporal_model1, name="model-test", ) dataset_fp.models = model print(dataset_fp) .. rst-class:: sphx-glr-script-out .. code-block:: none FluxPointsDataset ----------------- Name : dataset_lc Number of total flux points : 30 Number of fit bins : 25 Fit statistic type : chi2 Fit statistic value (-2 log(L)) : 1417.89 Number of models : 1 Number of parameters : 5 Number of free parameters : 3 Component 0: SkyModel Name : model-test Datasets names : None Spectral model type : PowerLawSpectralModel Spatial model type : Temporal model type : ExpDecayTemporalModel Parameters: index : 2.000 +/- 0.00 amplitude : 1.00e-12 +/- 0.0e+00 1 / (cm2 s TeV) reference (frozen): 1.000 TeV t0 : 10.000 +/- 0.00 h t_ref (frozen): 58909.000 d .. GENERATED FROM PYTHON SOURCE LINES 300-301 Fit the dataset .. GENERATED FROM PYTHON SOURCE LINES 301-306 .. code-block:: Python fit = Fit() result = fit.run(dataset_fp) 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 2.9966e+00 2.818e-02 ... False False amplitude 9.2654e-12 cm-2 s-1 TeV-1 3.260e-13 ... False True reference 1.0000e+00 TeV 0.000e+00 ... True False t0 6.2842e+00 h 2.538e-01 ... False False t_ref 5.8909e+04 d 0.000e+00 ... True False .. GENERATED FROM PYTHON SOURCE LINES 307-310 Now let’s plot model and data. We plot only the normalisation of the temporal model in relative units for one particular energy range .. GENERATED FROM PYTHON SOURCE LINES 310-313 .. code-block:: Python dataset_fp.plot_spectrum(axis_name="time") .. image-sg:: /tutorials/analysis-time/images/sphx_glr_light_curve_simulation_002.png :alt: light curve simulation :srcset: /tutorials/analysis-time/images/sphx_glr_light_curve_simulation_002.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 314-324 Fit the datasets ~~~~~~~~~~~~~~~~ Here, we apply the models directly to the simulated datasets. For modelling and fitting more complex flares, you should attach the relevant model to each group of ``datasets``. The parameters of a model in a given group of dataset will be tied. For more details on joint fitting in Gammapy, see the :doc:`/tutorials/analysis-3d/analysis_3d`. .. GENERATED FROM PYTHON SOURCE LINES 324-341 .. code-block:: Python # Define the model: spectral_model2 = PowerLawSpectralModel( index=2.0, amplitude="1e-12 cm-2 s-1 TeV-1", reference="1 TeV" ) temporal_model2 = ExpDecayTemporalModel(t0="10 h", t_ref=gti_t0.mjd * u.d) model2 = SkyModel( spectral_model=spectral_model2, temporal_model=temporal_model2, name="model-test2", ) display(model2.parameters.to_table()) datasets.models = model2 .. rst-class:: sphx-glr-script-out .. code-block:: none type name value unit error ... frozen is_norm link prior ---- --------- ---------- -------------- --------- ... ------ ------- ---- ----- index 2.0000e+00 0.000e+00 ... False False amplitude 1.0000e-12 cm-2 s-1 TeV-1 0.000e+00 ... False True reference 1.0000e+00 TeV 0.000e+00 ... True False t0 1.0000e+01 h 0.000e+00 ... False False t_ref 5.8909e+04 d 0.000e+00 ... True False .. GENERATED FROM PYTHON SOURCE LINES 342-343 Perform a joint fit .. GENERATED FROM PYTHON SOURCE LINES 343-348 .. code-block:: Python fit = Fit() result = fit.run(datasets=datasets) 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.0113e+00 3.271e-02 ... False False amplitude 9.4119e-12 cm-2 s-1 TeV-1 3.337e-13 ... False True reference 1.0000e+00 TeV 0.000e+00 ... True False t0 6.1101e+00 h 2.152e-01 ... False False t_ref 5.8909e+04 d 0.000e+00 ... True False .. GENERATED FROM PYTHON SOURCE LINES 349-352 We see that the fitted parameters are consistent between fitting flux points and datasets, and match well with the simulated ones .. GENERATED FROM PYTHON SOURCE LINES 355-364 Exercises --------- 1. Re-do the analysis with `~gammapy.datasets.MapDataset` instead of a `~gammapy.datasets.SpectrumDataset` 2. Model the flare of PKS 2155-304 which you obtained using the :doc:`/tutorials/analysis-time/light_curve_flare` tutorial. Use a combination of a Gaussian and Exponential flare profiles. 3. Do a joint fitting of the datasets. .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 18.651 seconds) .. _sphx_glr_download_tutorials_analysis-time_light_curve_simulation.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-time/light_curve_simulation.ipynb :alt: Launch binder :width: 150 px .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: light_curve_simulation.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: light_curve_simulation.py ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_