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

Multi instrument joint 3D and 1D analysis#

Prerequisites#

Context#

Some science studies require to combine heterogeneous data from various instruments to extract physical information. In particular, it is often useful to add flux measurements of a source at different energies to an analysis to better constrain the wide-band spectral parameters. This can be done using a joint fit of heterogeneous datasets.

Objectives: Constrain the spectral parameters of the gamma-ray emission from the Crab nebula between 10 GeV and 100 TeV, using a 3D Fermi dataset, a H.E.S.S. reduced spectrum and HAWC flux points.

Proposed approach#

This tutorial illustrates how to perform a joint modeling and fitting of the Crab Nebula spectrum using different datasets. The spectral parameters are optimized by combining a 3D analysis of Fermi-LAT data, a ON/OFF spectral analysis of HESS data, and flux points from HAWC.

In this tutorial we are going to use pre-made datasets. We prepared maps of the Crab region as seen by Fermi-LAT using the same event selection than the 3FHL catalog (7 years of data with energy from 10 GeV to 2 TeV). For the HESS ON/OFF analysis we used two observations from the first public data release with a significant signal from energy of about 600 GeV to 10 TeV. These observations have an offset of 0.5° and a zenith angle of 45-48°. The HAWC flux points data are taken from a recent analysis based on 2.5 years of data with energy between 300 Gev and 300 TeV.

The setup#

[1]:
from astropy import units as u
import matplotlib.pyplot as plt
from gammapy.modeling import Fit
from gammapy.modeling.models import Models, create_crab_spectral_model
from gammapy.datasets import Datasets, FluxPointsDataset, SpectrumDatasetOnOff
from gammapy.estimators import FluxPoints, FluxPointsEstimator
from gammapy.maps import MapAxis
from pathlib import Path

Data and models files#

The datasets serialization produce YAML files listing the datasets and models. In the following cells we show an example containning only the Fermi-LAT dataset and the Crab model.

Fermi-LAT-3FHL_datasets.yaml:

[2]:
!cat $GAMMAPY_DATA/fermi-3fhl-crab/Fermi-LAT-3FHL_datasets.yaml
datasets:
- name: Fermi-LAT
  type: MapDataset
  filename: Fermi-LAT-3FHL_data_Fermi-LAT.fits

We used as model a point source with a log-parabola spectrum. The initial parameters were taken from the latest Fermi-LAT catalog 4FGL, then we have re-optimized the spectral parameters for our dataset in the 10 GeV - 2 TeV energy range (fixing the source position).

Fermi-LAT-3FHL_models.yaml:

[3]:
!cat $GAMMAPY_DATA/fermi-3fhl-crab/Fermi-LAT-3FHL_models.yaml
components:
-   name: Crab Nebula
    type: SkyModel
    spectral:
        type: LogParabolaSpectralModel
        parameters:
        -   name: amplitude
            value: 0.018182745349064267
            unit: cm-2 s-1 TeV-1
            min: .nan
            max: .nan
            frozen: false
            error: 0.003026327991562108
        -   name: reference
            value: 5.054833602905273e-05
            unit: TeV
            min: .nan
            max: .nan
            frozen: true
            error: 0.0
        -   name: alpha
            value: 1.652368617859867
            unit: ''
            min: .nan
            max: .nan
            frozen: false
            error: 0.05762513693893088
        -   name: beta
            value: 0.03921700077803329
            unit: ''
            min: .nan
            max: .nan
            frozen: false
            error: 0.00521472221220211
    spatial:
        type: PointSpatialModel
        frame: icrs
        parameters:
        -   name: lon_0
            value: 83.63310241699219
            unit: deg
            min: .nan
            max: .nan
            frozen: true
            error: 0.0
        -   name: lat_0
            value: 22.019899368286133
            unit: deg
            min: -90.0
            max: 90.0
            frozen: true
            error: 0.0
-   type: FoVBackgroundModel
    datasets_names:
    - Fermi-LAT
    spectral:
        type: PowerLawNormSpectralModel
        parameters:
        -   name: norm
            value: 1.3004625872247901
            unit: ''
            min: 0.0
            max: .nan
            frozen: false
            error: 0.07512322002655547
        -   name: tilt
            value: 0.0
            unit: ''
            min: .nan
            max: .nan
            frozen: true
            error: 0.0
        -   name: reference
            value: 1.0
            unit: TeV
            min: .nan
            max: .nan
            frozen: true
            error: 0.0

Reading different datasets#

Fermi-LAT 3FHL: map dataset for 3D analysis#

For now we let’s use the datasets serialization only to read the 3D MapDataset associated to Fermi-LAT 3FHL data and models.

[4]:
path = Path("$GAMMAPY_DATA/fermi-3fhl-crab")
filename = path / "Fermi-LAT-3FHL_datasets.yaml"

datasets = Datasets.read(filename=filename)
[5]:
models = Models.read(path / "Fermi-LAT-3FHL_models.yaml")
print(models)
Models

Component 0: SkyModel

  Name                      : Crab Nebula
  Datasets names            : None
  Spectral model type       : LogParabolaSpectralModel
  Spatial  model type       : PointSpatialModel
  Temporal model type       :
  Parameters:
    amplitude                     :   1.82e-02   +/- 3.0e-03 1 / (cm2 s TeV)
    reference             (frozen):      0.000       TeV
    alpha                         :      1.652   +/-    0.06
    beta                          :      0.039   +/-    0.01
    lon_0                 (frozen):     83.633       deg
    lat_0                 (frozen):     22.020       deg

Component 1: FoVBackgroundModel

  Name                      : Fermi-LAT-bkg
  Datasets names            : ['Fermi-LAT']
  Spectral model type       : PowerLawNormSpectralModel
  Parameters:
    norm                          :      1.300   +/-    0.08
    tilt                  (frozen):      0.000
    reference             (frozen):      1.000       TeV


We get the Crab model in order to share it with the other datasets

[6]:
print(models["Crab Nebula"])
SkyModel

  Name                      : Crab Nebula
  Datasets names            : None
  Spectral model type       : LogParabolaSpectralModel
  Spatial  model type       : PointSpatialModel
  Temporal model type       :
  Parameters:
    amplitude                     :   1.82e-02   +/- 3.0e-03 1 / (cm2 s TeV)
    reference             (frozen):      0.000       TeV
    alpha                         :      1.652   +/-    0.06
    beta                          :      0.039   +/-    0.01
    lon_0                 (frozen):     83.633       deg
    lat_0                 (frozen):     22.020       deg


HESS-DL3: 1D ON/OFF dataset for spectral fitting#

The ON/OFF datasets can be read from PHA files following the OGIP standards. We read the PHA files from each observation, and compute a stacked dataset for simplicity. Then the Crab spectral model previously defined is added to the dataset.

[7]:
datasets_hess = Datasets()

for obs_id in [23523, 23526]:
    dataset = SpectrumDatasetOnOff.read(
        f"$GAMMAPY_DATA/joint-crab/spectra/hess/pha_obs{obs_id}.fits"
    )
    datasets_hess.append(dataset)

dataset_hess = datasets_hess.stack_reduce(name="HESS")

datasets.append(dataset_hess)
[8]:
print(datasets)
Datasets
--------

Dataset 0:

  Type       : MapDataset
  Name       : Fermi-LAT
  Instrument :
  Models     :

Dataset 1:

  Type       : SpectrumDatasetOnOff
  Name       : HESS
  Instrument :
  Models     :


HAWC: 1D dataset for flux point fitting#

The HAWC flux point are taken from https://arxiv.org/pdf/1905.12518.pdf Then these flux points are read from a pre-made FITS file and passed to a FluxPointsDataset together with the source spectral model.

[9]:
# read flux points from https://arxiv.org/pdf/1905.12518.pdf
filename = "$GAMMAPY_DATA/hawc_crab/HAWC19_flux_points.fits"
flux_points_hawc = FluxPoints.read(
    filename, reference_model=create_crab_spectral_model("meyer")
)

dataset_hawc = FluxPointsDataset(data=flux_points_hawc, name="HAWC")

datasets.append(dataset_hawc)
[10]:
print(datasets)
Datasets
--------

Dataset 0:

  Type       : MapDataset
  Name       : Fermi-LAT
  Instrument :
  Models     :

Dataset 1:

  Type       : SpectrumDatasetOnOff
  Name       : HESS
  Instrument :
  Models     :

Dataset 2:

  Type       : FluxPointsDataset
  Name       : HAWC
  Instrument :
  Models     :


Datasets serialization#

The datasets object contains each dataset previously defined. It can be saved on disk as datasets.yaml, models.yaml, and several data files specific to each dataset. Then the datasets can be rebuild later from these files.

[11]:
path = Path("crab-3datasets")
path.mkdir(exist_ok=True)

filename = path / "crab_10GeV_100TeV_datasets.yaml"

datasets.write(filename, overwrite=True)
HDU 'MASK_FIT' not found
[12]:
!cat crab-3datasets/crab_10GeV_100TeV_datasets.yaml
datasets:
- name: Fermi-LAT
  type: MapDataset
  filename: Fermi-LAT.fits
- name: HESS
  type: SpectrumDatasetOnOff
  filename: pha_obsHESS.fits
- name: HAWC
  type: FluxPointsDataset
  filename: HAWC.fits
[13]:
datasets = Datasets.read(filename)
datasets.models = models
No reference model set for FluxMaps. Assuming point source with E^-2 spectrum.
[14]:
print(datasets)
Datasets
--------

Dataset 0:

  Type       : MapDataset
  Name       : Fermi-LAT
  Instrument :
  Models     : ['Crab Nebula', 'Fermi-LAT-bkg']

Dataset 1:

  Type       : SpectrumDatasetOnOff
  Name       : HESS
  Instrument :
  Models     : ['Crab Nebula']

Dataset 2:

  Type       : FluxPointsDataset
  Name       : HAWC
  Instrument :
  Models     : ['Crab Nebula']


Joint analysis#

We run the fit on the Datasets object that include a dataset for each instrument

[15]:
%%time
fit_joint = Fit()
results_joint = fit_joint.run(datasets=datasets)
print(results_joint)
HDU 'MASK_FIT' not found
OptimizeResult

        backend    : minuit
        method     : migrad
        success    : True
        message    : Optimization terminated successfully.
        nfev       : 302
        total stat : -12697.22

CovarianceResult

        backend    : minuit
        method     : hesse
        success    : True
        message    : Hesse terminated successfully.

CPU times: user 2.25 s, sys: 47.4 ms, total: 2.3 s
Wall time: 2.74 s

Let’s display only the parameters of the Crab spectral model

[16]:
crab_spec = datasets[0].models["Crab Nebula"].spectral_model
print(crab_spec)
LogParabolaSpectralModel

  type      name     value         unit        error   min max frozen is_norm link
-------- --------- ---------- -------------- --------- --- --- ------ ------- ----
spectral amplitude 3.9740e-03 cm-2 s-1 TeV-1 3.125e-04 nan nan  False    True
spectral reference 5.0548e-05            TeV 0.000e+00 nan nan   True   False
spectral     alpha 1.2634e+00                1.707e-02 nan nan  False   False
spectral      beta 6.1321e-02                9.448e-04 nan nan  False   False

We can compute flux points for Fermi-LAT and HESS datasets in order plot them together with the HAWC flux point.

[17]:
# compute Fermi-LAT and HESS flux points
energy_edges = MapAxis.from_energy_bounds("10 GeV", "2 TeV", nbin=5).edges

flux_points_fermi = FluxPointsEstimator(
    energy_edges=energy_edges,
    source="Crab Nebula",
).run([datasets["Fermi-LAT"]])


energy_edges = MapAxis.from_bounds(
    1, 15, nbin=6, interp="log", unit="TeV"
).edges

flux_points_hess = FluxPointsEstimator(
    energy_edges=energy_edges, source="Crab Nebula", selection_optional=["ul"]
).run([datasets["HESS"]])

Now, Let’s plot the Crab spectrum fitted and the flux points of each instrument.

[18]:
# display spectrum and flux points
plt.figure(figsize=(8, 6))

energy_bounds = [0.01, 300] * u.TeV
sed_type = "e2dnde"

ax = crab_spec.plot(
    energy_bounds=energy_bounds, sed_type=sed_type, label="Model"
)
crab_spec.plot_error(ax=ax, energy_bounds=energy_bounds, sed_type=sed_type)

flux_points_fermi.plot(ax=ax, sed_type=sed_type, label="Fermi-LAT")
flux_points_hess.plot(ax=ax, sed_type=sed_type, label="HESS")
flux_points_hawc.plot(ax=ax, sed_type=sed_type, label="HAWC")

ax.set_xlim(energy_bounds)
plt.legend();
../../../_images/tutorials_analysis_3D_analysis_mwl_30_0.png
[ ]: