.. 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-105 .. code-block:: Python from pathlib import Path # Check package versions 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 106-111 Setup ----- As usual, we’ll start with some setup … .. GENERATED FROM PYTHON SOURCE LINES 111-134 .. 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 135-137 Check setup ----------- .. GENERATED FROM PYTHON SOURCE LINES 137-143 .. 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 : /Users/mregeard/Workspace/dev/code/gammapy/gammapy/.tox/build_docs/bin/python python_version : 3.11.10 machine : x86_64 system : Darwin Gammapy package: version : 1.3.dev1205+g00f44f94ac path : /Users/mregeard/Workspace/dev/code/gammapy/gammapy/.tox/build_docs/lib/python3.11/site-packages/gammapy Other packages: numpy : 1.26.4 scipy : 1.14.1 astropy : 5.2.2 regions : 0.10 click : 8.1.7 yaml : 6.0.2 IPython : 8.28.0 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.2 ray : 2.37.0 Gammapy environment variables: GAMMAPY_DATA : /Users/mregeard/Workspace/dev/code/gammapy/gammapy-data/ .. GENERATED FROM PYTHON SOURCE LINES 144-152 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 152-158 .. 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 159-166 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 166-172 .. 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 173-185 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 185-201 .. 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 202-207 Run data reduction chain ------------------------ We begin with the configuration of the maker classes: .. GENERATED FROM PYTHON SOURCE LINES 207-224 .. 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 225-235 .. 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 236-239 Plot off regions ---------------- .. GENERATED FROM PYTHON SOURCE LINES 239-247 .. 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 /Users/mregeard/Workspace/dev/code/gammapy/gammapy/.tox/build_docs/lib/python3.11/site-packages/regions/shapes/circle.py:160: UserWarning: Setting the 'color' property will override the edgecolor or facecolor properties. return Circle(xy=xy, radius=radius, **mpl_kwargs) /Users/mregeard/Workspace/dev/code/gammapy/gammapy/.tox/build_docs/lib/python3.11/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 248-254 Source statistic ---------------- Next we’re going to look at the overall source statistics in our signal region. .. GENERATED FROM PYTHON SOURCE LINES 254-259 .. 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 sqrt_ts background ... stat_sum counts_off acceptance acceptance_off alpha ... ------- ------ ----------------- ------------------ ------------------ ... ------------------ ---------- ---------- ------------------ ------------------- stacked 149 139.25 20.449683569684254 9.75 ... 433.5372460592368 117 18.0 216.0 0.0833333358168602 stacked 303 280.75 28.446255766130005 22.250001907348633 ... 823.6909484036685 267 37.0 443.99993896484375 0.0833333432674408 stacked 439 408.7743835449219 36.17588855633591 30.225610733032227 ... 1325.2636829165463 594 56.0 1100.523681640625 0.05088486522436142 stacked 550 512.135498046875 40.90086300803346 37.864498138427734 ... 1701.2360501021876 869 74.0 1698.3189697265625 0.04357249662280083 .. GENERATED FROM PYTHON SOURCE LINES 260-261 And make the corresponding plots .. GENERATED FROM PYTHON SOURCE LINES 261-287 .. 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 288-293 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 293-301 .. 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 302-304 If you want to read back the datasets from disk you can use: .. GENERATED FROM PYTHON SOURCE LINES 304-312 .. 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 313-321 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 321-339 .. 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 340-343 Fit quality and model residuals ------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 346-348 We can access the results dictionary to see if the fit converged: .. GENERATED FROM PYTHON SOURCE LINES 348-352 .. 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 353-355 and check the best-fit parameters .. GENERATED FROM PYTHON SOURCE LINES 355-359 .. code-block:: Python display(result_joint.models.to_parameters_table()) .. rst-class:: sphx-glr-script-out .. code-block:: none model type name value unit error min max frozen link prior ----- ---- --------- ---------- -------------- --------- --- --- ------ ---- ----- crab index 2.2727e+00 1.566e-01 nan nan False crab amplitude 4.7913e-11 cm-2 s-1 TeV-1 3.600e-12 nan nan False crab reference 1.0000e+00 TeV 0.000e+00 nan nan True crab lambda_ 1.2097e-01 TeV-1 5.382e-02 nan nan False crab alpha 1.0000e+00 0.000e+00 nan nan True .. GENERATED FROM PYTHON SOURCE LINES 360-363 A simple way to inspect the model residuals is using the function `~SpectrumDataset.plot_fit()` .. GENERATED FROM PYTHON SOURCE LINES 363-369 .. code-block:: Python ax_spectrum, ax_residuals = datasets[0].plot_fit() ax_spectrum.set_ylim(0.1, 40) 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 370-373 For more ways of assessing fit quality, please refer to the dedicated :doc:`/tutorials/api/fitting` tutorial. .. GENERATED FROM PYTHON SOURCE LINES 376-385 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 create an instance of the `~gammapy.estimators.FluxPointsEstimator`, by passing the dataset and the energy binning: .. GENERATED FROM PYTHON SOURCE LINES 385-392 .. code-block:: Python fpe = FluxPointsEstimator( energy_edges=energy_axis.edges, source="crab", selection_optional="all" ) flux_points = fpe.run(datasets=datasets) .. GENERATED FROM PYTHON SOURCE LINES 393-395 Here is a the table of the resulting flux points: .. GENERATED FROM PYTHON SOURCE LINES 395-399 .. 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 dnde dnde_err dnde_errp dnde_errn ... stat_null stat_scan is_ul counts success norm_scan TeV TeV TeV 1 / (cm2 s TeV) 1 / (cm2 s TeV) 1 / (cm2 s TeV) 1 / (cm2 s TeV) ... ------ ------ ------ --------------- --------------- --------------- --------------- ... --------- ------------------ ----- ------------ ------- -------------- 0.112 0.100 0.125 nan nan nan nan ... 0.000 nan .. nan False 0.0 .. 0.0 False 0.200 .. 5.000 0.139 0.125 0.156 nan nan nan nan ... 0.000 nan .. nan False 0.0 .. 0.0 False 0.200 .. 5.000 0.174 0.156 0.195 nan nan nan nan ... 0.000 nan .. nan False 0.0 .. 0.0 False 0.200 .. 5.000 0.217 0.195 0.243 nan nan nan nan ... 0.000 nan .. nan False 0.0 .. 0.0 False 0.200 .. 5.000 0.271 0.243 0.303 nan nan nan nan ... 0.000 nan .. nan False 0.0 .. 0.0 False 0.200 .. 5.000 0.339 0.303 0.379 nan nan nan nan ... 0.000 nan .. nan False 0.0 .. 0.0 False 0.200 .. 5.000 0.423 0.379 0.473 nan nan nan nan ... 0.000 nan .. nan False 0.0 .. 0.0 False 0.200 .. 5.000 0.528 0.473 0.590 nan nan nan nan ... 0.000 nan .. nan False 0.0 .. 0.0 False 0.200 .. 5.000 0.659 0.590 0.737 1.132e-10 1.632e-11 1.702e-11 1.564e-11 ... 140.184 62.683 .. 252.535 False 0.0 .. 0.0 True 0.200 .. 5.000 0.823 0.737 0.920 6.950e-11 7.224e-12 7.462e-12 7.000e-12 ... 326.560 140.739 .. 451.124 False 30.0 .. 25.0 True 0.200 .. 5.000 1.028 0.920 1.148 3.822e-11 4.560e-12 4.732e-12 4.391e-12 ... 260.522 105.409 .. 381.859 False 30.0 .. 17.0 True 0.200 .. 5.000 1.283 1.148 1.434 2.164e-11 2.930e-12 3.055e-12 2.808e-12 ... 201.897 77.194 .. 311.795 False 13.0 .. 18.0 True 0.200 .. 5.000 1.602 1.434 1.790 1.527e-11 2.062e-12 2.150e-12 1.975e-12 ... 210.704 88.719 .. 221.364 False 19.0 .. 14.0 True 0.200 .. 5.000 2.000 1.790 2.235 7.852e-12 1.237e-12 1.299e-12 1.178e-12 ... 145.510 62.080 .. 203.491 False 18.0 .. 7.0 True 0.200 .. 5.000 2.497 2.235 2.790 4.494e-12 7.803e-13 8.236e-13 7.386e-13 ... 130.346 57.228 .. 169.536 False 4.0 .. 7.0 True 0.200 .. 5.000 3.117 2.790 3.483 2.674e-12 5.110e-13 5.419e-13 4.811e-13 ... 101.683 45.318 .. 123.910 False 9.0 .. 10.0 True 0.200 .. 5.000 3.892 3.483 4.348 8.756e-13 2.576e-13 2.800e-13 2.362e-13 ... 40.269 13.594 .. 125.092 False 8.0 .. 2.0 True 0.200 .. 5.000 4.859 4.348 5.429 7.435e-13 1.952e-13 2.118e-13 1.795e-13 ... 58.151 31.556 .. 80.689 False 10.0 .. 3.0 True 0.200 .. 5.000 6.066 5.429 6.778 3.779e-13 1.197e-13 1.320e-13 1.081e-13 ... 40.590 15.545 .. 51.420 False 2.0 .. 3.0 True 0.200 .. 5.000 7.573 6.778 8.462 2.240e-13 8.066e-14 8.979e-14 7.205e-14 ... 28.471 16.394 .. 34.290 False 4.0 .. 1.0 True 0.200 .. 5.000 9.454 8.462 10.564 1.190e-13 5.279e-14 5.981e-14 4.659e-14 ... 18.617 12.081 .. 21.493 False 1.0 .. 4.0 True 0.200 .. 5.000 11.803 10.564 13.189 5.214e-14 3.008e-14 3.529e-14 2.546e-14 ... 15.845 12.285 .. 18.679 False 0.0 .. 0.0 True 0.200 .. 5.000 14.736 13.189 16.465 7.642e-15 1.124e-14 1.505e-14 nan ... 4.027 3.401 .. 13.778 True 0.0 .. 0.0 True 0.200 .. 5.000 18.397 16.465 20.556 1.494e-14 1.127e-14 1.407e-14 8.796e-15 ... 7.853 6.221 .. 3.778 False 1.0 .. 0.0 True 0.200 .. 5.000 22.968 20.556 25.663 6.571e-19 5.174e-17 2.837e-15 nan ... 0.254 0.423 .. 4.486 True 0.0 .. 0.0 True 0.200 .. 5.000 28.675 25.663 32.040 -1.615e-31 1.798e-17 1.892e-15 nan ... 0.366 0.443 .. 2.288 True 0.0 .. 0.0 True 0.200 .. 5.000 35.799 32.040 40.000 2.836e-21 1.502e-18 nan nan ... 0.110 0.139 .. 0.855 True 0.0 .. 0.0 True 0.200 .. 5.000 .. GENERATED FROM PYTHON SOURCE LINES 400-403 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 403-410 .. 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") ax.set_xlim(0.6, 40) 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 411-415 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 418-421 The final plot with the best fit model, flux points and residuals can be quickly made like this: .. GENERATED FROM PYTHON SOURCE LINES 421-430 .. code-block:: Python flux_points_dataset = FluxPointsDataset( data=flux_points, models=model_best_joint.copy() ) ax, _ = flux_points_dataset.plot_fit() ax.set_xlim(0.6, 40) 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 431-438 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 438-442 .. code-block:: Python dataset_stacked = Datasets(datasets).stack_reduce() .. GENERATED FROM PYTHON SOURCE LINES 443-447 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 447-457 .. 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 458-459 And display the parameter table .. GENERATED FROM PYTHON SOURCE LINES 459-465 .. 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 min max frozen link prior ---- --------- ---------- -------------- --------- --- --- ------ ---- ----- index 2.2727e+00 1.566e-01 nan nan False amplitude 4.7913e-11 cm-2 s-1 TeV-1 3.600e-12 nan nan False reference 1.0000e+00 TeV 0.000e+00 nan nan True lambda_ 1.2097e-01 TeV-1 5.382e-02 nan nan False alpha 1.0000e+00 0.000e+00 nan nan True type name value unit error min max frozen link prior ---- --------- ---------- -------------- --------- --- --- ------ ---- ----- index 2.2785e+00 1.563e-01 nan nan False amplitude 4.7800e-11 cm-2 s-1 TeV-1 3.566e-12 nan nan False reference 1.0000e+00 TeV 0.000e+00 nan nan True lambda_ 1.1830e-01 TeV-1 5.329e-02 nan nan False alpha 1.0000e+00 0.000e+00 nan nan True .. GENERATED FROM PYTHON SOURCE LINES 466-470 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 470-498 .. 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 499-515 A note on statistics -------------------- Different statistic are available for the FluxPointDataset : * chi2 : estimate from chi2 statistics. * profile : estimate from interpolation of the likelihood profile. * distrib : estimate from probability distributions, assuming that flux points correspond to asymmetric gaussians and upper limits complementary error functions. Default is `chi2`, in that case upper limits are ignored and the mean of asymetrics error is used. So it is recommended to use `profile` if `stat_scan` is available on flux points. The `distrib` case provides an approximation if the `profile` is not available which allows to take into accounts upper limit and asymetrics error. In the example below we can see that the `profile` case matches exactly the result from the joint analysis of the ON/OFF datasets using `wstat` (as labelled). .. GENERATED FROM PYTHON SOURCE LINES 515-552 .. code-block:: Python def plot_stat(fp_dataset): fig, ax = plt.subplots() plot_kwargs = { "energy_bounds": [0.1, 30] * u.TeV, "sed_type": "e2dnde", "ax": ax, } fp_dataset.data.plot(energy_power=2, ax=ax) model_best_joint.spectral_model.plot( color="b", lw=0.5, **plot_kwargs, label="wstat" ) stat_types = ["chi2", "profile", "distrib"] colors = ["red", "g", "c"] lss = ["--", ":", "--"] for ks, stat in enumerate(stat_types): fp_dataset.stat_type = stat fit = Fit() fit.run([fp_dataset]) fp_dataset.models[0].spectral_model.plot( color=colors[ks], ls=lss[ks], **plot_kwargs, label=stat ) fp_dataset.models[0].spectral_model.plot_error( facecolor=colors[ks], **plot_kwargs ) plt.legend() plot_stat(flux_points_dataset) .. image-sg:: /tutorials/analysis-1d/images/sphx_glr_spectral_analysis_008.png :alt: spectral analysis :srcset: /tutorials/analysis-1d/images/sphx_glr_spectral_analysis_008.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 553-573 .. code-block:: Python # In order to avoid discrepancies due to the treatment of upper limits # we can utilise the `~gammapy.estimators.utils.resample_energy_edges` # for defining energy bins in which the minimum number of `sqrt_ts` is 2. # In that case all the statistics definitions give equivalent results. # energy_edges = resample_energy_edges(dataset_stacked, conditions={"sqrt_ts_min": 2}) fpe_no_ul = FluxPointsEstimator( energy_edges=energy_edges, source="crab", selection_optional="all" ) flux_points_no_ul = fpe_no_ul.run(datasets=datasets) flux_points_dataset_no_ul = FluxPointsDataset( data=flux_points_no_ul, models=model_best_joint.copy(), ) plot_stat(flux_points_dataset_no_ul) .. image-sg:: /tutorials/analysis-1d/images/sphx_glr_spectral_analysis_009.png :alt: spectral analysis :srcset: /tutorials/analysis-1d/images/sphx_glr_spectral_analysis_009.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 574-589 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 592-603 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 22.172 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 ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: spectral_analysis.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_