.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "tutorials/analysis-1d/spectral_analysis.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-1d_spectral_analysis.py: Spectral analysis ================= Perform a full region based on-off spectral analysis and fit the resulting datasets. Prerequisites ------------- - Understanding how spectral extraction is performed in Cherenkov astronomy, in particular regarding OFF background measurements. - Understanding the basics data reduction and modeling/fitting process with the gammapy library API as shown in the :doc:`/tutorials/starting/analysis_2` tutorial. Context ------- While 3D analyses allow in principle to consider complex field of views containing overlapping gamma-ray sources, in many cases we might have an observation with a single, strong, point-like source in the field of view. A spectral analysis, in that case, might consider all the events inside a source (or ON) region and bin them in energy only, obtaining 1D datasets. In classical Cherenkov astronomy, the background estimation technique associated with this method measures the number of events in OFF regions taken in regions of the field-of-view devoid of gamma-ray emitters, where the background rate is assumed to be equal to the one in the ON region. This allows to use a specific fit statistics for ON-OFF measurements, the wstat (see `~gammapy.stats.wstat`), where no background model is assumed. Background is treated as a set of nuisance parameters. This removes some systematic effects connected to the choice or the quality of the background model. But this comes at the expense of larger statistical uncertainties on the fitted model parameters. **Objective: perform a full region based spectral analysis of 4 Crab observations of H.E.S.S. data release 1 and fit the resulting datasets.** Introduction ------------ Here, as usual, we use the `~gammapy.data.DataStore` to retrieve a list of selected observations (`~gammapy.data.Observations`). Then, we define the ON region containing the source and the geometry of the `~gammapy.datasets.SpectrumDataset` object we want to produce. We then create the corresponding dataset Maker. We have to define the Maker object that will extract the OFF counts from reflected regions in the field-of-view. To ensure we use data in an energy range where the quality of the IRFs is good enough we also create a safe range Maker. We can then proceed with data reduction with a loop over all selected observations to produce datasets in the relevant geometry. We can then explore the resulting datasets and look at the cumulative signal and significance of our source. We finally proceed with model fitting. In practice, we have to: - Create a `~gammapy.data.DataStore` pointing to the relevant data - Apply an observation selection to produce a list of observations, a `~gammapy.data.Observations` object. - Define a geometry of the spectrum we want to produce: - Create a `~regions.CircleSkyRegion` for the ON extraction region - Create a `~gammapy.maps.MapAxis` for the energy binnings: one for the reconstructed (i.e.measured) energy, the other for the true energy (i.e.the one used by IRFs and models) - Create the necessary makers : - the spectrum dataset maker : `~gammapy.makers.SpectrumDatasetMaker` - the OFF background maker, here a `~gammapy.makers.ReflectedRegionsBackgroundMaker` - and the safe range maker : `~gammapy.makers.SafeMaskMaker` - Perform the data reduction loop. And for every observation: - Apply the makers sequentially to produce a `~gammapy.datasets.SpectrumDatasetOnOff` - Append it to list of datasets - Define the `~gammapy.modeling.models.SkyModel` to apply to the dataset. - Create a `~gammapy.modeling.Fit` object and run it to fit the model parameters - Apply a `~gammapy.estimators.FluxPointsEstimator` to compute flux points for the spectral part of the fit. .. GENERATED FROM PYTHON SOURCE LINES 94-106 .. code-block:: Python from pathlib import Path # Check package versions import numpy as np import astropy.units as u from astropy.coordinates import Angle, SkyCoord from regions import CircleSkyRegion # %matplotlib inline import matplotlib.pyplot as plt .. GENERATED FROM PYTHON SOURCE LINES 107-112 Setup ----- As usual, we’ll start with some setup … .. GENERATED FROM PYTHON SOURCE LINES 112-135 .. code-block:: Python from IPython.display import display from gammapy.data import DataStore from gammapy.datasets import ( Datasets, FluxPointsDataset, SpectrumDataset, SpectrumDatasetOnOff, ) from gammapy.estimators import FluxPointsEstimator from gammapy.estimators.utils import resample_energy_edges from gammapy.makers import ( ReflectedRegionsBackgroundMaker, SafeMaskMaker, SpectrumDatasetMaker, ) from gammapy.maps import MapAxis, RegionGeom, WcsGeom from gammapy.modeling import Fit from gammapy.modeling.models import ( ExpCutoffPowerLawSpectralModel, SkyModel, create_crab_spectral_model, ) .. GENERATED FROM PYTHON SOURCE LINES 136-138 Check setup ----------- .. GENERATED FROM PYTHON SOURCE LINES 138-144 .. code-block:: Python from gammapy.utils.check import check_tutorials_setup from gammapy.visualization import plot_spectrum_datasets_off_regions 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.dev296+g18b65d5cc 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.8.4 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.20.0 Gammapy environment variables: GAMMAPY_DATA : /home/runner/work/gammapy-docs/gammapy-docs/gammapy-datasets/dev .. GENERATED FROM PYTHON SOURCE LINES 145-153 Load Data --------- First, we select and load some H.E.S.S. observations of the Crab nebula. We will access the events, effective area, energy dispersion, livetime and PSF for containment correction. .. GENERATED FROM PYTHON SOURCE LINES 153-159 .. code-block:: Python datastore = DataStore.from_dir("$GAMMAPY_DATA/hess-dl3-dr1/") obs_ids = [23523, 23526, 23559, 23592] observations = datastore.get_observations(obs_ids) .. GENERATED FROM PYTHON SOURCE LINES 160-167 Define Target Region -------------------- The next step is to define a signal extraction region, also known as on region. In the simplest case this is just a `CircleSkyRegion `__. .. GENERATED FROM PYTHON SOURCE LINES 167-173 .. code-block:: Python target_position = SkyCoord(ra=83.63, dec=22.01, unit="deg", frame="icrs") on_region_radius = Angle("0.11 deg") on_region = CircleSkyRegion(center=target_position, radius=on_region_radius) .. GENERATED FROM PYTHON SOURCE LINES 174-186 Create exclusion mask --------------------- We will use the reflected regions method to place off regions to estimate the background level in the on region. To make sure the off regions don’t contain gamma-ray emission, we create an exclusion mask. Using http://gamma-sky.net/ we find that there’s only one known gamma-ray source near the Crab nebula: the AGN called `RGB J0521+212 `__ at GLON = 183.604 deg and GLAT = -8.708 deg. .. GENERATED FROM PYTHON SOURCE LINES 186-202 .. code-block:: Python exclusion_region = CircleSkyRegion( center=SkyCoord(183.604, -8.708, unit="deg", frame="galactic"), radius=0.5 * u.deg, ) skydir = target_position.galactic geom = WcsGeom.create( npix=(150, 150), binsz=0.05, skydir=skydir, proj="TAN", frame="icrs" ) exclusion_mask = ~geom.region_mask([exclusion_region]) exclusion_mask.plot() plt.show() .. image-sg:: /tutorials/analysis-1d/images/sphx_glr_spectral_analysis_001.png :alt: spectral analysis :srcset: /tutorials/analysis-1d/images/sphx_glr_spectral_analysis_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 203-208 Run data reduction chain ------------------------ We begin with the configuration of the maker classes: .. GENERATED FROM PYTHON SOURCE LINES 208-225 .. code-block:: Python energy_axis = MapAxis.from_energy_bounds( 0.1, 40, nbin=10, per_decade=True, unit="TeV", name="energy" ) energy_axis_true = MapAxis.from_energy_bounds( 0.05, 100, nbin=20, per_decade=True, unit="TeV", name="energy_true" ) geom = RegionGeom.create(region=on_region, axes=[energy_axis]) dataset_empty = SpectrumDataset.create(geom=geom, energy_axis_true=energy_axis_true) dataset_maker = SpectrumDatasetMaker( containment_correction=True, selection=["counts", "exposure", "edisp"] ) bkg_maker = ReflectedRegionsBackgroundMaker(exclusion_mask=exclusion_mask) safe_mask_masker = SafeMaskMaker(methods=["aeff-max"], aeff_percent=10) .. GENERATED FROM PYTHON SOURCE LINES 226-236 .. code-block:: Python datasets = Datasets() for obs_id, observation in zip(obs_ids, observations): dataset = dataset_maker.run(dataset_empty.copy(name=str(obs_id)), observation) dataset_on_off = bkg_maker.run(dataset, observation) dataset_on_off = safe_mask_masker.run(dataset_on_off, observation) datasets.append(dataset_on_off) print(datasets) .. rst-class:: sphx-glr-script-out .. code-block:: none Datasets -------- Dataset 0: Type : SpectrumDatasetOnOff Name : 23523 Instrument : HESS Models : Dataset 1: Type : SpectrumDatasetOnOff Name : 23526 Instrument : HESS Models : Dataset 2: Type : SpectrumDatasetOnOff Name : 23559 Instrument : HESS Models : Dataset 3: Type : SpectrumDatasetOnOff Name : 23592 Instrument : HESS Models : .. GENERATED FROM PYTHON SOURCE LINES 237-240 Plot off regions ---------------- .. GENERATED FROM PYTHON SOURCE LINES 240-248 .. code-block:: Python plt.figure() ax = exclusion_mask.plot() on_region.to_pixel(ax.wcs).plot(ax=ax, edgecolor="k") plot_spectrum_datasets_off_regions(ax=ax, datasets=datasets) plt.show() .. image-sg:: /tutorials/analysis-1d/images/sphx_glr_spectral_analysis_002.png :alt: spectral analysis :srcset: /tutorials/analysis-1d/images/sphx_glr_spectral_analysis_002.png :class: sphx-glr-single-img .. 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/regions/shapes/circle.py:161: UserWarning: Setting the 'color' property will override the edgecolor or facecolor properties. return Circle(xy=xy, radius=radius, **mpl_kwargs) /home/runner/work/gammapy-docs/gammapy-docs/gammapy/.tox/build_docs/lib/python3.9/site-packages/gammapy/visualization/datasets.py:84: UserWarning: Setting the 'color' property will override the edgecolor or facecolor properties. handle = Patch(**plot_kwargs) .. GENERATED FROM PYTHON SOURCE LINES 249-255 Source statistic ---------------- Next we’re going to look at the overall source statistics in our signal region. .. GENERATED FROM PYTHON SOURCE LINES 255-260 .. code-block:: Python info_table = datasets.info_table(cumulative=True) display(info_table) .. rst-class:: sphx-glr-script-out .. code-block:: none name counts excess ... acceptance_off alpha ... ------- ------ ----------------- ... ------------------ ------------------- stacked 149 139.25 ... 216.0 0.0833333358168602 stacked 303 280.75 ... 227.99996948242188 0.0833333432674408 stacked 439 408.7743835449219 ... 373.3919677734375 0.05088486522436142 stacked 550 512.135498046875 ... 436.05487060546875 0.04357249662280083 .. GENERATED FROM PYTHON SOURCE LINES 261-262 And make the corresponding plots .. GENERATED FROM PYTHON SOURCE LINES 262-288 .. code-block:: Python fig, (ax_excess, ax_sqrt_ts) = plt.subplots(figsize=(10, 4), ncols=2, nrows=1) ax_excess.plot( info_table["livetime"].to("h"), info_table["excess"], marker="o", ls="none", ) ax_excess.set_title("Excess") ax_excess.set_xlabel("Livetime [h]") ax_excess.set_ylabel("Excess events") ax_sqrt_ts.plot( info_table["livetime"].to("h"), info_table["sqrt_ts"], marker="o", ls="none", ) ax_sqrt_ts.set_title("Sqrt(TS)") ax_sqrt_ts.set_xlabel("Livetime [h]") ax_sqrt_ts.set_ylabel("Sqrt(TS)") plt.show() .. image-sg:: /tutorials/analysis-1d/images/sphx_glr_spectral_analysis_003.png :alt: Excess, Sqrt(TS) :srcset: /tutorials/analysis-1d/images/sphx_glr_spectral_analysis_003.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 289-294 Finally you can write the extracted datasets to disk using the OGIP format (PHA, ARF, RMF, BKG, see `here `__ for details): .. GENERATED FROM PYTHON SOURCE LINES 294-302 .. code-block:: Python path = Path("spectrum_analysis") path.mkdir(exist_ok=True) for dataset in datasets: dataset.write(filename=path / f"obs_{dataset.name}.fits.gz", overwrite=True) .. GENERATED FROM PYTHON SOURCE LINES 303-305 If you want to read back the datasets from disk you can use: .. GENERATED FROM PYTHON SOURCE LINES 305-313 .. code-block:: Python datasets = Datasets() for obs_id in obs_ids: filename = path / f"obs_{obs_id}.fits.gz" datasets.append(SpectrumDatasetOnOff.read(filename)) .. GENERATED FROM PYTHON SOURCE LINES 314-322 Fit spectrum ------------ Now we’ll fit a global model to the spectrum. First we do a joint likelihood fit to all observations. If you want to stack the observations see below. We will also produce a debug plot in order to show how the global fit matches one of the individual observations. .. GENERATED FROM PYTHON SOURCE LINES 322-340 .. code-block:: Python spectral_model = ExpCutoffPowerLawSpectralModel( amplitude=1e-12 * u.Unit("cm-2 s-1 TeV-1"), index=2, lambda_=0.1 * u.Unit("TeV-1"), reference=1 * u.TeV, ) model = SkyModel(spectral_model=spectral_model, name="crab") datasets.models = [model] fit_joint = Fit() result_joint = fit_joint.run(datasets=datasets) # we make a copy here to compare it later model_best_joint = model.copy() .. GENERATED FROM PYTHON SOURCE LINES 341-344 Fit quality and model residuals ------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 347-349 We can access the results dictionary to see if the fit converged: .. GENERATED FROM PYTHON SOURCE LINES 349-353 .. 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 : 244 total stat : 86.12 CovarianceResult backend : minuit method : hesse success : True message : Hesse terminated successfully. .. GENERATED FROM PYTHON SOURCE LINES 354-356 and check the best-fit parameters .. GENERATED FROM PYTHON SOURCE LINES 356-360 .. code-block:: Python display(result_joint.models.to_parameters_table()) .. rst-class:: sphx-glr-script-out .. code-block:: none model type name value unit ... max frozen is_norm link prior ----- ---- --------- ---------- -------------- ... --- ------ ------- ---- ----- crab index 2.2727e+00 ... nan False False crab amplitude 4.7913e-11 cm-2 s-1 TeV-1 ... nan False True crab reference 1.0000e+00 TeV ... nan True False crab lambda_ 1.2097e-01 TeV-1 ... nan False False crab alpha 1.0000e+00 ... nan True False .. GENERATED FROM PYTHON SOURCE LINES 361-364 A simple way to inspect the model residuals is using the function `~SpectrumDataset.plot_fit()` .. GENERATED FROM PYTHON SOURCE LINES 364-371 .. code-block:: Python ax_spectrum, ax_residuals = datasets[0].plot_fit() ax_spectrum.set_ylim(0.1, 40) datasets[0].plot_masks(ax=ax_spectrum) plt.show() .. image-sg:: /tutorials/analysis-1d/images/sphx_glr_spectral_analysis_004.png :alt: spectral analysis :srcset: /tutorials/analysis-1d/images/sphx_glr_spectral_analysis_004.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 372-375 For more ways of assessing fit quality, please refer to the dedicated :doc:`/tutorials/api/fitting` tutorial. .. GENERATED FROM PYTHON SOURCE LINES 378-387 Compute Flux Points ------------------- To round up our analysis we can compute flux points by fitting the norm of the global model in energy bands. We can utilise the `~gammapy.estimators.utils.resample_energy_edges` for defining the energy bins in which the minimum number of `sqrt_ts` is 2. To do so we first stack the individual datasets, only for obtaining the energies: .. GENERATED FROM PYTHON SOURCE LINES 387-392 .. code-block:: Python dataset_stacked = Datasets(datasets).stack_reduce() energy_edges = resample_energy_edges(dataset_stacked, conditions={"sqrt_ts_min": 2}) .. GENERATED FROM PYTHON SOURCE LINES 393-397 Now we create an instance of the `~gammapy.estimators.FluxPointsEstimator`, by passing the dataset and the energy binning: .. GENERATED FROM PYTHON SOURCE LINES 397-404 .. code-block:: Python fpe = FluxPointsEstimator( energy_edges=energy_edges, source="crab", selection_optional="all" ) flux_points = fpe.run(datasets=datasets) .. GENERATED FROM PYTHON SOURCE LINES 405-407 Here is a the table of the resulting flux points: .. GENERATED FROM PYTHON SOURCE LINES 407-411 .. code-block:: Python display(flux_points.to_table(sed_type="dnde", formatted=True)) .. rst-class:: sphx-glr-script-out .. code-block:: none e_ref e_min e_max ... success norm_scan stat_scan TeV TeV TeV ... ------ ------ ------ ... ------- -------------- ------------------ 0.659 0.590 0.737 ... True 0.200 .. 5.000 62.683 .. 252.535 0.823 0.737 0.920 ... True 0.200 .. 5.000 140.739 .. 451.124 1.028 0.920 1.148 ... True 0.200 .. 5.000 105.409 .. 381.859 1.283 1.148 1.434 ... True 0.200 .. 5.000 77.194 .. 311.795 1.602 1.434 1.790 ... True 0.200 .. 5.000 88.719 .. 221.364 2.000 1.790 2.235 ... True 0.200 .. 5.000 62.080 .. 203.491 2.497 2.235 2.790 ... True 0.200 .. 5.000 57.228 .. 169.536 3.117 2.790 3.483 ... True 0.200 .. 5.000 45.318 .. 123.910 3.892 3.483 4.348 ... True 0.200 .. 5.000 13.594 .. 125.092 4.859 4.348 5.429 ... True 0.200 .. 5.000 31.556 .. 80.689 6.066 5.429 6.778 ... True 0.200 .. 5.000 15.545 .. 51.420 7.573 6.778 8.462 ... True 0.200 .. 5.000 16.394 .. 34.290 9.454 8.462 10.564 ... True 0.200 .. 5.000 12.081 .. 21.493 20.556 10.564 40.000 ... True 0.200 .. 5.000 22.913 .. 43.864 .. GENERATED FROM PYTHON SOURCE LINES 412-415 Now we plot the flux points and their likelihood profiles. For the plotting of upper limits we choose a threshold of TS < 4. .. GENERATED FROM PYTHON SOURCE LINES 415-421 .. code-block:: Python fig, ax = plt.subplots() flux_points.plot(ax=ax, sed_type="e2dnde", color="darkorange") flux_points.plot_ts_profiles(ax=ax, sed_type="e2dnde") plt.show() .. image-sg:: /tutorials/analysis-1d/images/sphx_glr_spectral_analysis_005.png :alt: spectral analysis :srcset: /tutorials/analysis-1d/images/sphx_glr_spectral_analysis_005.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 422-426 Note: it is also possible to plot the flux distribution with the spectral model overlaid, but you must ensure the axis binning is identical for the flux points and integral flux. .. GENERATED FROM PYTHON SOURCE LINES 429-432 The final plot with the best fit model, flux points and residuals can be quickly made like this: .. GENERATED FROM PYTHON SOURCE LINES 432-438 .. code-block:: Python flux_points_dataset = FluxPointsDataset(data=flux_points, models=model_best_joint) flux_points_dataset.plot_fit() plt.show() .. image-sg:: /tutorials/analysis-1d/images/sphx_glr_spectral_analysis_006.png :alt: spectral analysis :srcset: /tutorials/analysis-1d/images/sphx_glr_spectral_analysis_006.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 439-446 Stack observations ------------------ An alternative approach to fitting the spectrum is stacking all observations first and the fitting a model. For this we first stack the individual datasets: .. GENERATED FROM PYTHON SOURCE LINES 446-450 .. code-block:: Python dataset_stacked = Datasets(datasets).stack_reduce() .. GENERATED FROM PYTHON SOURCE LINES 451-455 Again we set the model on the dataset we would like to fit (in this case it’s only a single one) and pass it to the `~gammapy.modeling.Fit` object: .. GENERATED FROM PYTHON SOURCE LINES 455-465 .. code-block:: Python dataset_stacked.models = model stacked_fit = Fit() result_stacked = stacked_fit.run([dataset_stacked]) # Make a copy to compare later model_best_stacked = model.copy() print(result_stacked) .. rst-class:: sphx-glr-script-out .. code-block:: none OptimizeResult backend : minuit method : migrad success : True message : Optimization terminated successfully. nfev : 54 total stat : 8.16 CovarianceResult backend : minuit method : hesse success : True message : Hesse terminated successfully. .. GENERATED FROM PYTHON SOURCE LINES 466-467 And display the parameter table .. GENERATED FROM PYTHON SOURCE LINES 467-473 .. code-block:: Python display(model_best_joint.parameters.to_table()) display(model_best_stacked.parameters.to_table()) .. rst-class:: sphx-glr-script-out .. code-block:: none type name value unit error ... frozen is_norm link prior ---- --------- ---------- -------------- --------- ... ------ ------- ---- ----- index 2.2727e+00 1.566e-01 ... False False amplitude 4.7913e-11 cm-2 s-1 TeV-1 3.600e-12 ... False True reference 1.0000e+00 TeV 0.000e+00 ... True False lambda_ 1.2097e-01 TeV-1 5.382e-02 ... False False alpha 1.0000e+00 0.000e+00 ... True False type name value unit error ... frozen is_norm link prior ---- --------- ---------- -------------- --------- ... ------ ------- ---- ----- index 2.2785e+00 1.563e-01 ... False False amplitude 4.7800e-11 cm-2 s-1 TeV-1 3.566e-12 ... False True reference 1.0000e+00 TeV 0.000e+00 ... True False lambda_ 1.1830e-01 TeV-1 5.329e-02 ... False False alpha 1.0000e+00 0.000e+00 ... True False .. GENERATED FROM PYTHON SOURCE LINES 474-478 Finally, we compare the results of our stacked analysis to a previously published Crab Nebula Spectrum for reference. This is available in `~gammapy.modeling.models.create_crab_spectral_model`. .. GENERATED FROM PYTHON SOURCE LINES 478-506 .. code-block:: Python fig, ax = plt.subplots() plot_kwargs = { "energy_bounds": [0.1, 30] * u.TeV, "sed_type": "e2dnde", "yunits": u.Unit("erg cm-2 s-1"), "ax": ax, } # plot stacked model model_best_stacked.spectral_model.plot(**plot_kwargs, label="Stacked analysis result") model_best_stacked.spectral_model.plot_error(facecolor="blue", alpha=0.3, **plot_kwargs) # plot joint model model_best_joint.spectral_model.plot( **plot_kwargs, label="Joint analysis result", ls="--" ) model_best_joint.spectral_model.plot_error(facecolor="orange", alpha=0.3, **plot_kwargs) create_crab_spectral_model("hess_ecpl").plot( **plot_kwargs, label="Crab reference", ) ax.legend() plt.show() # sphinx_gallery_thumbnail_number = 5 .. image-sg:: /tutorials/analysis-1d/images/sphx_glr_spectral_analysis_007.png :alt: spectral analysis :srcset: /tutorials/analysis-1d/images/sphx_glr_spectral_analysis_007.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 507-522 Exercises --------- Now you have learned the basics of a spectral analysis with Gammapy. To practice you can continue with the following exercises: - Fit a different spectral model to the data. You could try `~gammapy.modeling.models.ExpCutoffPowerLawSpectralModel` or `~gammapy.modeling.models.LogParabolaSpectralModel`. - Compute flux points for the stacked dataset. - Create a `~gammapy.datasets.FluxPointsDataset` with the flux points you have computed for the stacked dataset and fit the flux points again with obe of the spectral models. How does the result compare to the best fit model, that was directly fitted to the counts data? .. GENERATED FROM PYTHON SOURCE LINES 525-536 What next? ---------- The methods shown in this tutorial is valid for point-like or midly extended sources where we can assume that the IRF taken at the region center is valid over the whole region. If one wants to extract the 1D spectrum of a large source and properly average the response over the extraction region, one has to use a different approach explained in the :doc:`/tutorials/analysis-1d/extended_source_spectral_analysis` tutorial. .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 15.044 seconds) .. _sphx_glr_download_tutorials_analysis-1d_spectral_analysis.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-1d/spectral_analysis.ipynb :alt: Launch binder :width: 150 px .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: spectral_analysis.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: spectral_analysis.py ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_