.. 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.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 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 create an instance of the `~gammapy.estimators.FluxPointsEstimator`, by passing the dataset and the energy binning: .. GENERATED FROM PYTHON SOURCE LINES 387-394 .. 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 395-397 Here is a the table of the resulting flux points: .. GENERATED FROM PYTHON SOURCE LINES 397-401 .. 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.112 0.100 0.125 ... False 0.200 .. 5.000 nan .. nan 0.139 0.125 0.156 ... False 0.200 .. 5.000 nan .. nan 0.174 0.156 0.195 ... False 0.200 .. 5.000 nan .. nan 0.217 0.195 0.243 ... False 0.200 .. 5.000 nan .. nan 0.271 0.243 0.303 ... False 0.200 .. 5.000 nan .. nan 0.339 0.303 0.379 ... False 0.200 .. 5.000 nan .. nan 0.423 0.379 0.473 ... False 0.200 .. 5.000 nan .. nan 0.528 0.473 0.590 ... False 0.200 .. 5.000 nan .. nan 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 ... ... ... ... ... ... ... 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 11.803 10.564 13.189 ... True 0.200 .. 5.000 12.285 .. 18.679 14.736 13.189 16.465 ... True 0.200 .. 5.000 3.401 .. 13.778 18.397 16.465 20.556 ... True 0.200 .. 5.000 6.221 .. 3.778 22.968 20.556 25.663 ... False 0.200 .. 5.000 0.423 .. 4.486 28.675 25.663 32.040 ... True 0.200 .. 5.000 0.443 .. 2.288 35.799 32.040 40.000 ... False 0.200 .. 5.000 0.139 .. 0.855 Length = 27 rows .. GENERATED FROM PYTHON SOURCE LINES 402-405 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 405-412 .. 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 413-417 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 420-423 The final plot with the best fit model, flux points and residuals can be quickly made like this: .. GENERATED FROM PYTHON SOURCE LINES 423-432 .. 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 433-440 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 440-444 .. code-block:: Python dataset_stacked = Datasets(datasets).stack_reduce() .. GENERATED FROM PYTHON SOURCE LINES 445-449 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 449-459 .. 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 460-461 And display the parameter table .. GENERATED FROM PYTHON SOURCE LINES 461-467 .. 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 468-472 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 472-500 .. 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 501-517 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 517-555 .. 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 556-576 .. 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 577-592 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 595-606 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 29.131 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 `_