forked from open-gamma-ray-astro/joint-crab
-
Notifications
You must be signed in to change notification settings - Fork 0
/
conf.py
159 lines (137 loc) · 4.81 KB
/
conf.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
"""Configuration and data access helper class"""
import numpy as np
import astropy.units as u
from astropy.coordinates import SkyCoord
from astropy.table import Table
from gammapy.spectrum import SpectrumObservationList
from gammapy.maps import Map
DATASETS = [
{
"name": "fermi",
"label": "Fermi-LAT",
"obs_ids": [0],
"on_radius": "0.3 deg",
"containment_correction": True,
"energy_range": {"min": "0.03 TeV", "max": "2 TeV"},
"color": "#21ABCD",
},
{
"name": "magic",
"label": "MAGIC",
"on_radius": "0.141 deg",
"containment_correction": False,
"energy_range": {"min": "0.08 TeV", "max": "30 TeV"},
"color": "#FF9933",
},
{
"name": "veritas",
"label": "VERITAS",
"on_radius": "0.1 deg",
"containment_correction": False,
"energy_range": {"min": "0.15 TeV", "max": "30 TeV"},
"color": "#893F45",
},
{
"name": "fact",
"label": "FACT",
"on_radius": "0.1732 deg",
"containment_correction": False,
"energy_range": {"min": "0.4 TeV", "max": "30 TeV"},
"color": "#3EB489",
},
{
"name": "hess",
"label": "H.E.S.S.",
"on_radius": "0.11 deg",
"containment_correction": True,
"energy_range": {"min": "0.66 TeV", "max": "30 TeV"},
"color": "#002E63",
},
{
"name": "joint",
"label": "joint fit",
"energy_range": {"min": "0.03 TeV", "max": "30 TeV"},
"color": "crimson",
},
]
class PlotConfig:
fontsize = 15
fontsize_contours = 18
label_energy = r"$E\,/\,\mathrm{TeV}$"
e2dnde_label = (
r"$E^2 \cdot {\rm d}\phi/{\rm d}E\,/\,({\rm erg}\,{\rm cm}^{-2} {\rm s}^{-1})$"
)
class Config:
"""Configuration options."""
source_pos = SkyCoord("83d37m59.0988s", "22d00m52.2s")
plot = PlotConfig()
all_datasets = ["fermi", "magic", "veritas", "fact", "hess"]
all_datasets_plus_joint = all_datasets + ["joint"]
def __init__(self):
self.datasets = {_["name"]: Dataset.from_dict(_) for _ in DATASETS}
# energy bins for all the datasets 10 GeV to 100 TeV, 20 bins per decade
self.energy_bins = np.logspace(np.log10(0.01), np.log10(100), 81) * u.TeV
@staticmethod
def get_exclusion_mask():
return Map.read("results/maps/exclusion_mask.fits.gz")
class Dataset(object):
"""Information about each dataset."""
def __init__(
self,
name,
label,
obs_ids,
on_radius,
containment_correction,
energy_range,
color,
):
self.name = name
self.label = label
self.obs_ids = obs_ids
self.on_radius = on_radius
self.containment_correction = containment_correction
self.energy_range = energy_range
self.color = color
# directories to read / store results per dataset
self.data_path = f"data/{self.name}"
self.spectra_path = f"results/spectra/{self.name}"
@classmethod
def from_dict(cls, d):
"""function to initialize the class starting from a dictionary with keys
dictionary as defined per each source in the config.yaml
"""
e_min = u.Quantity(d["energy_range"]["min"])
e_max = u.Quantity(d["energy_range"]["max"])
energy_range = [e_min.value, e_max.value] * e_min.unit
if d["name"] in ["magic", "veritas", "fact", "hess"]:
table = Table.read(f"data/{d['name']}/obs-index.fits.gz")
obs_ids = np.array(table["OBS_ID"])
else:
obs_ids = None
return cls(
name=d["name"],
label=d["label"],
obs_ids=obs_ids,
on_radius=u.Quantity(d.get("on_radius", "nan")),
containment_correction=d.get("containment_correction", "nan"),
energy_range=energy_range,
color=d["color"],
)
def get_energy_mask(self, energy):
return (self.energy_range[0] <= energy) & (energy < self.energy_range[1])
def get_SpectrumObservationList(self):
"""load the OGIP files and return a SpectrumObservationList
SpectrumObservationList has already a method to read from a directory
http://docs.gammapy.org/dev/api/gammapy.spectrum.SpectrumObservationList.html
"""
if self.name == "joint":
spec_obs_list = SpectrumObservationList()
# extend the list adding all the other SpectrumObservationList
for name in {"fermi", "magic", "veritas", "fact", "hess"}:
spec_obs = SpectrumObservationList.read(f"results/spectra/{name}")
spec_obs_list.extend(spec_obs)
else:
spec_obs_list = SpectrumObservationList.read(self.spectra_path)
return spec_obs_list
config = Config()