Note
Click here to download the full example code or to run this example in your browser via Binder
Multi instrument joint 3D and 1D analysis#
Joint 3D analysis using 3D Fermi datasets, a H.E.S.S. reduced spectrum and HAWC flux points.
Prerequisites#
Handling of Fermi-LAT data with Gammapy see the Fermi-LAT with Gammapy tutorial.
Knowledge of spectral analysis to produce 1D On-Off datasets, see the following doc:
/tutorials/analysis-1d/spectral_analysistutorial.Using flux points to directly fit a model (without forward-folding) from the doc:
/tutorials/analysis-1d/sed_fittingtutorial.
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#
from pathlib import Path
from astropy import units as u
import matplotlib.pyplot as plt
from gammapy.datasets import Datasets, FluxPointsDataset, SpectrumDatasetOnOff
from gammapy.estimators import FluxPoints, FluxPointsEstimator
from gammapy.maps import MapAxis
from gammapy.modeling import Fit
from gammapy.modeling.models import Models, create_crab_spectral_model
Check setup#
from gammapy.utils.check import check_tutorials_setup
from gammapy.utils.scripts import make_path
check_tutorials_setup()
System:
python_executable : /Users/terrier/Code/anaconda3/envs/gammapy-dev/bin/python
python_version : 3.8.13
machine : x86_64
system : Darwin
Gammapy package:
version : 1.0rc2
path : /Users/terrier/Code/gammapy-dev/gammapy/gammapy
Other packages:
numpy : 1.22.4
scipy : 1.9.3
astropy : 5.1
regions : 0.6
click : 8.1.3
yaml : 6.0
IPython : 8.4.0
jupyterlab : 3.4.8
matplotlib : 3.5.3
pandas : 1.5.0
healpy : 1.16.1
iminuit : 2.17.0
sherpa : 4.15.0
naima : 0.10.0
emcee : 3.1.3
corner : 2.2.1
Gammapy environment variables:
GAMMAPY_DATA : /Users/terrier/Code/gammapy-dev/gammapy-data
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:
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:
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.
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
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.
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)
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.
# 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)
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.
path = Path("crab-3datasets")
path.mkdir(exist_ok=True)
filename = path / "crab_10GeV_100TeV_datasets.yaml"
datasets.write(filename, overwrite=True)
datasets = Datasets.read(filename)
datasets.models = models
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
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.
Let’s display only the parameters of the Crab spectral model
LogParabolaSpectralModel
type name value unit ... max frozen is_norm link
-------- --------- ---------- -------------- ... --- ------ ------- ----
spectral amplitude 3.9740e-03 cm-2 s-1 TeV-1 ... nan False True
spectral reference 5.0548e-05 TeV ... nan True False
spectral alpha 1.2634e+00 ... nan False False
spectral beta 6.1321e-02 ... nan False False
We can compute flux points for Fermi-LAT and HESS datasets in order plot them together with the HAWC flux point.
# 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.
# display spectrum and flux points
fig, ax = plt.subplots(figsize=(8, 6))
energy_bounds = [0.01, 300] * u.TeV
sed_type = "e2dnde"
crab_spec.plot(ax=ax, 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)
ax.legend()
plt.show()
