.. 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 : /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 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_maker = 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_maker.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 /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 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 ... acceptance_off alpha ... ------- ------ ----------------- ... ------------------ ------------------- stacked 149 139.25 ... 216.0 0.0833333358168602 stacked 303 280.75 ... 443.99993896484375 0.0833333432674408 stacked 439 408.7743835449219 ... 1100.523681640625 0.05088486522436142 stacked 550 512.135498046875 ... 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 ... min max frozen link prior ----- ---- --------- ---------- -------------- ... --- --- ------ ---- ----- crab index 2.2727e+00 ... nan nan False crab amplitude 4.7913e-11 cm-2 s-1 TeV-1 ... nan nan False crab reference 1.0000e+00 TeV ... nan nan True crab lambda_ 1.2097e-01 TeV-1 ... nan nan False crab alpha 1.0000e+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 ... counts success norm_scan TeV TeV TeV 1 / (cm2 s TeV) ... ------ ------ ------ --------------- ... ------------ ------- -------------- 0.112 0.100 0.125 nan ... 0.0 .. 0.0 False 0.200 .. 5.000 0.139 0.125 0.156 nan ... 0.0 .. 0.0 False 0.200 .. 5.000 0.174 0.156 0.195 nan ... 0.0 .. 0.0 False 0.200 .. 5.000 0.217 0.195 0.243 nan ... 0.0 .. 0.0 False 0.200 .. 5.000 0.271 0.243 0.303 nan ... 0.0 .. 0.0 False 0.200 .. 5.000 0.339 0.303 0.379 nan ... 0.0 .. 0.0 False 0.200 .. 5.000 0.423 0.379 0.473 nan ... 0.0 .. 0.0 False 0.200 .. 5.000 0.528 0.473 0.590 nan ... 0.0 .. 0.0 False 0.200 .. 5.000 0.659 0.590 0.737 1.132e-10 ... 0.0 .. 0.0 True 0.200 .. 5.000 0.823 0.737 0.920 6.950e-11 ... 30.0 .. 25.0 True 0.200 .. 5.000 ... ... ... ... ... ... ... ... 4.859 4.348 5.429 7.435e-13 ... 10.0 .. 3.0 True 0.200 .. 5.000 6.066 5.429 6.778 3.779e-13 ... 2.0 .. 3.0 True 0.200 .. 5.000 7.573 6.778 8.462 2.240e-13 ... 4.0 .. 1.0 True 0.200 .. 5.000 9.454 8.462 10.564 1.190e-13 ... 1.0 .. 4.0 True 0.200 .. 5.000 11.803 10.564 13.189 5.214e-14 ... 0.0 .. 0.0 True 0.200 .. 5.000 14.736 13.189 16.465 7.642e-15 ... 0.0 .. 0.0 True 0.200 .. 5.000 18.397 16.465 20.556 1.494e-14 ... 1.0 .. 0.0 True 0.200 .. 5.000 22.968 20.556 25.663 -5.332e-31 ... 0.0 .. 0.0 True 0.200 .. 5.000 28.675 25.663 32.040 7.267e-23 ... 0.0 .. 0.0 True 0.200 .. 5.000 35.799 32.040 40.000 -4.119e-32 ... 0.0 .. 0.0 True 0.200 .. 5.000 Length = 27 rows .. 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-517 A note on statistics -------------------- Different statistic are available for the `~gammapy.datasets.FluxPointsDataset` : - 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 517-554 .. 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 555-560 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. .. GENERATED FROM PYTHON SOURCE LINES 561-575 .. code-block:: Python 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 576-591 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 594-605 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 31.521 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/v1.3?urlpath=lab/tree/notebooks/1.3/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 `_