This is a fixed-text formatted version of a Jupyter notebook

Flux Profile Estimation

This tutorial shows how to estimate flux profiles.

Prerequisites

Knowledge of 3D data reduction and datasets used in Gammapy, see for instance the first analysis tutorial.

Context

A useful tool to study and compare the saptial distribution of flux in images and data cubes is the measurement of flxu profiles. Flux profiles can show spatial correlations of gamma-ray data with e.g. gas maps or other type of gamma-ray data. Most commonly flux profiles are measured along some preferred coordinate axis, either radially distance from a source of interest, along longitude and latitude coordinate axes or along the path defined by two spatial coordinates.

Proposed Approach

Flux profile estimation essentially works by estimating flux points for a set of predefined spatially connected regions. For radial flux profiles the shape of the regions are annuli with a common center, for linear profiles it’s typically a rectangular shape.

We will work on a pre-computed MapDataset of Fermi-LAT data, use Region to define the structure of the bins of the flux profile and run the actually profile extraction using the FluxProfileEstimator

Setup and Imports

[1]:
%matplotlib inline
import matplotlib.pyplot as plt
from astropy import units as u
from astropy.coordinates import SkyCoord
import numpy as np
[2]:
from gammapy.datasets import MapDataset
from gammapy.estimators import FluxProfileEstimator, FluxPoints
from gammapy.maps import RegionGeom
from gammapy.utils.regions import (
    make_concentric_annulus_sky_regions,
    make_orthogonal_rectangle_sky_regions,
)
from gammapy.modeling.models import PowerLawSpectralModel

Read and Introduce Data

[3]:
dataset = MapDataset.read(
    "$GAMMAPY_DATA/fermi-3fhl-gc/fermi-3fhl-gc.fits.gz", name="fermi-dataset"
)

This is what the counts image we will work with looks like:

[4]:
counts_image = dataset.counts.sum_over_axes()
counts_image.smooth("0.1 deg").plot(stretch="sqrt");
../../../_images/tutorials_analysis_3D_flux_profiles_8_0.png

There are 400x200 pixels in the dataset and 11 energy bins between 10 GeV and 2 TeV:

[5]:
print(dataset.counts)
WcsNDMap

        geom  : WcsGeom
        axes  : ['lon', 'lat', 'energy']
        shape : (400, 200, 11)
        ndim  : 3
        unit  :
        dtype : >i4

Profile Estimation

Configuration

We start by defining a list of spatially connected regions along the galactic longitude axis. For this there is a helper function make_orthogonal_rectangle_sky_regions. The individual region bins for the profile have a height of 3 deg and in total there are 31 bins. The starts from lon = 10 deg tand goes to lon = 350 deg. In addition we have to specify the wcs to take into account possible projections effects on the region definition:

[6]:
regions = make_orthogonal_rectangle_sky_regions(
    start_pos=SkyCoord("10d", "0d", frame="galactic"),
    end_pos=SkyCoord("350d", "0d", frame="galactic"),
    wcs=counts_image.geom.wcs,
    height="3 deg",
    nbin=51,
)

We can use the RegionGeom object to illustrate the regions on top of the counts image:

[7]:
geom = RegionGeom.create(region=regions)
ax = counts_image.smooth("0.1 deg").plot(stretch="sqrt")
geom.plot_region(ax=ax, color="w");
../../../_images/tutorials_analysis_3D_flux_profiles_14_0.png

Next we create the FluxProfileEstimator. For the estimation of the flux profile we assume a spectral model with a power-law shape and an index of 2.3

[8]:
flux_profile_estimator = FluxProfileEstimator(
    regions=regions,
    spectrum=PowerLawSpectralModel(index=2.3),
    energy_edges=[10, 2000] * u.GeV,
    selection_optional=["ul"],
)

We can see the full configuration by printing the estimator object:

[9]:
print(flux_profile_estimator)
FluxProfileEstimator
--------------------

  energy_edges           : [  10. 2000.] GeV
  fit                    : <gammapy.modeling.fit.Fit object at 0x7f274464e250>
  n_sigma                : 1
  n_sigma_ul             : 2
  norm_max               : 5
  norm_min               : 0.2
  norm_n_values          : 11
  norm_values            : None
  null_value             : 0
  reoptimize             : False
  selection_optional     : ['ul']
  source                 : 0
  spectrum               : PowerLawSpectralModel
  sum_over_energy_groups : False

Run Estimation

Now we can run the profile estimation and explore the results:

[10]:
%%time
profile = flux_profile_estimator.run(datasets=dataset)
CPU times: user 21.9 s, sys: 50.9 ms, total: 22 s
Wall time: 22 s
[11]:
print(profile)
FluxPoints
----------

  geom                   : RegionGeom
  axes                   : ['lon', 'lat', 'energy', 'projected-distance']
  shape                  : (1, 1, 1, 51)
  quantities             : ['norm', 'norm_err', 'norm_ul', 'ts', 'npred', 'npred_excess', 'stat', 'counts', 'success']
  ref. model             : pl
  n_sigma                : 1
  n_sigma_ul             : 2
  sqrt_ts_threshold_ul   : 2
  sed type init          : likelihood

We can see the flux profile is represented by a FluxPoints object with a projected-distance axis, which defines the main axis the flux profile is measured along. The lon and lat axes can be ignored.

Plotting Results

Let us directly plot the result using FluxPoints.plot():

[12]:
ax = profile.plot(sed_type="dnde")
ax.set_yscale("linear")
../../../_images/tutorials_analysis_3D_flux_profiles_23_0.png

Based on the spectral model we specified above we can also plot in any other sed type, e.g. energy flux and define a different threshold when to plot upper limits:

[13]:
profile.sqrt_ts_threshold_ul = 2

ax = profile.plot(sed_type="eflux")
ax.set_yscale("linear")
../../../_images/tutorials_analysis_3D_flux_profiles_25_0.png

We can also plot any other quantity of interest, that is defined on the FluxPoints result object. E.g. the predicted total counts, background counts and excess counts:

[14]:
quantities = ["npred", "npred_excess", "npred_background"]

for quantity in quantities:
    profile[quantity].plot(label=quantity.title())

plt.ylabel("Counts ")
[14]:
Text(0, 0.5, 'Counts ')
../../../_images/tutorials_analysis_3D_flux_profiles_27_1.png

Serialisation and I/O

The profile can be serialised using FluxPoints.write(), given a specific format:

[15]:
profile.write(
    filename="flux_profile_fermi.fits",
    format="profile",
    overwrite=True,
    sed_type="dnde",
)
[16]:
profile_new = FluxPoints.read(
    filename="flux_profile_fermi.fits", format="profile"
)
No reference model set for FluxMaps. Assuming point source with E^-2 spectrum.
[17]:
ax = profile_new.plot()
ax.set_yscale("linear")
../../../_images/tutorials_analysis_3D_flux_profiles_31_0.png

The profile can be serialised to a ~astropy.table.Table object using:

[18]:
table = profile.to_table(format="profile", formatted=True)
table
[18]:
Table length=51
x_minx_maxx_refe_ref [1]e_min [1]e_max [1]ref_dnde [1]ref_flux [1]ref_eflux [1]norm [1]norm_err [1]norm_ul [1]ts [1]sqrt_ts [1]npred [1,1]npred_excess [1,1]stat [1]is_ul [1]counts [1,1]success [1]
degdegdegkeVkeVkeV1 / (cm2 s TeV)keV / (cm2 s TeV)keV2 / (cm2 s TeV)
float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float32float64boolfloat64bool
-0.19607843137254920.19607843137254920.070710678.11910000000.000500000000.0004.428e-103.043e-019.166e+06nannannannannannan0.00.000False0.0False
0.19607843137254920.58823529411764670.39215686274509870710678.11910000000.000500000000.0004.428e-103.043e-019.166e+06nannannannannannan0.00.000False0.0False
0.58823529411764670.98039215686274430.784313725490195570710678.11910000000.000500000000.0004.428e-103.043e-019.166e+060.1720.1240.4342.0591.435162.2304898612695517.448376-818.848True163.0True
0.98039215686274431.37254901960784361.17647058823529470710678.11910000000.000500000000.0004.428e-103.043e-019.166e+061.8940.2072.322120.59110.981450.20181889164763192.60909-3013.446False448.0True
1.37254901960784361.7647058823529421.568627450980392870710678.11910000000.000500000000.0004.428e-103.043e-019.166e+063.1460.2403.639280.64816.753601.6822045018936320.15994-4348.966False599.0True
1.7647058823529422.156862745098041.96078431372549170710678.11910000000.000500000000.0004.428e-103.043e-019.166e+060.7360.1831.11618.9914.358357.665030487863274.95968-2218.623False354.0True
2.156862745098042.5490196078431382.352941176470588870710678.11910000000.000500000000.0004.428e-103.043e-019.166e+060.6560.1861.04214.2143.770363.833381365701366.89659-2452.717False367.0True
2.5490196078431382.9411764705882362.74509803921568770710678.11910000000.000500000000.0004.428e-103.043e-019.166e+060.4500.1790.8217.0112.648339.4818273315226545.943058-2174.005False339.0True
2.9411764705882363.33333333333333443.13725490196078570710678.11910000000.000500000000.0004.428e-103.043e-019.166e+062.4210.2252.885174.14613.196528.0599574498646247.276-3887.531False531.0True
............................................................
15.88235294117646616.27450980392159616.0784313725490370710678.11910000000.000500000000.0004.428e-103.043e-019.166e+060.4270.1700.7806.9792.642324.690197080841944.554733-1957.576False321.0True
16.27450980392159616.66666666666669616.47058823529414470710678.11910000000.000500000000.0004.428e-103.043e-019.166e+060.7600.1951.16317.5534.190428.0052065967902379.45124-2767.621False421.0True
16.66666666666669617.05882352941177516.86274509803923470710678.11910000000.000500000000.0004.428e-103.043e-019.166e+060.6750.1971.08213.2773.644437.1341506314762570.588745-2830.574False431.0True
17.05882352941177517.450980392156917.25490196078433870710678.11910000000.000500000000.0004.428e-103.043e-019.166e+060.2560.1830.6342.0571.434381.2475096218942726.727816-2356.744True374.0True
17.450980392156917.84313725490200417.64705882352945270710678.11910000000.000500000000.0004.428e-103.043e-019.166e+060.1740.1820.5510.9510.975366.4533491178318.242867-2463.364True370.0True
17.84313725490200418.23529411764708318.03921568627454570710678.11910000000.000500000000.0004.428e-103.043e-019.166e+060.5930.1920.99010.7053.272408.055049644950662.12828-2803.387False410.0True
18.23529411764708318.62745098039215818.4313725490196270710678.11910000000.000500000000.0004.428e-103.043e-019.166e+060.3060.1730.6663.3421.828334.14141424656932.074398-2186.572True336.0True
18.62745098039215819.0196078431372618.8235294117647170710678.11910000000.000500000000.0004.428e-103.043e-019.166e+060.0100.1240.2690.0050.071171.811564914664730.9959075-877.032True172.0True
19.0196078431372619.4117647058823919.21568627450982470710678.11910000000.000500000000.0004.428e-103.043e-019.166e+06nannannannannannan0.00.000False0.0False
19.4117647058823919.80392156862749419.60784313725494270710678.11910000000.000500000000.0004.428e-103.043e-019.166e+06nannannannannannan0.00.000False0.0False

No we can also estimate a radial profile starting from the Galactic center:

[19]:
regions = make_concentric_annulus_sky_regions(
    center=SkyCoord("0d", "0d", frame="galactic"),
    radius_max="1.5 deg",
    nbin=11,
)

Again we first illustrate the regions:

[20]:
geom = RegionGeom.create(region=regions)
gc_image = counts_image.cutout(
    position=SkyCoord("0d", "0d", frame="galactic"), width=3 * u.deg
)
ax = gc_image.smooth("0.1 deg").plot(stretch="sqrt")
geom.plot_region(ax=ax, color="w");
/usr/share/miniconda/envs/gammapy-dev/lib/python3.7/site-packages/scipy/optimize/_numdiff.py:579: RuntimeWarning: invalid value encountered in true_divide
  J_transposed[i] = df / dx
../../../_images/tutorials_analysis_3D_flux_profiles_37_1.png

This time we define two energy bins and include the fit statistic profile in the computation:

[21]:
flux_profile_estimator = FluxProfileEstimator(
    regions=regions,
    spectrum=PowerLawSpectralModel(index=2.3),
    energy_edges=[10, 100, 2000] * u.GeV,
    selection_optional=["ul", "scan"],
    norm_values=np.linspace(-1, 5, 11),
)
[22]:
profile = flux_profile_estimator.run(datasets=dataset)

We can directly plot the result:

[23]:
ax = profile.plot(axis_name="projected-distance", sed_type="flux")
../../../_images/tutorials_analysis_3D_flux_profiles_42_0.png

However because of the powerlaw spectrum the flux at high energies is much lower. To extract the profile at high energies only we can use:

[24]:
profile_high = profile.slice_by_idx({"energy": slice(1, 2)})

And now plot the points together with the likelihood profiles:

[25]:
ax = profile_high.plot(sed_type="eflux", color="tab:orange")
profile_high.plot_ts_profiles(ax=ax, sed_type="eflux")
ax.set_yscale("linear")
../../../_images/tutorials_analysis_3D_flux_profiles_46_0.png
[ ]: