Mercurial > repos > astroteam > hess_astro_tool

comparison Spectrum_gammapy.py @ 1:593c4b45eda5 draft default tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit 2f23ec010eab8d33d2760f061286756e17015af7
author astroteam
date Thu, 18 Apr 2024 09:26:15 +0000
parents 02e4bb4fa10c
children
comparison
equal deleted inserted replaced
0:02e4bb4fa10c 1:593c4b45eda5
1 #!/usr/bin/env python
2 # coding: utf-8
3
4 # flake8: noqa
5
6 import json
7 import os
8 import shutil
9
10 import astropy.units as u
11 import matplotlib.pyplot as plt
12 import numpy as np
13 from astropy.coordinates import Angle, SkyCoord
14 from astropy.time import Time
15 from gammapy.data import DataStore
16 from gammapy.datasets import (
17 Datasets,
18 FluxPointsDataset,
19 SpectrumDataset,
20 SpectrumDatasetOnOff,
21 )
22 from gammapy.makers import (
23 ReflectedRegionsBackgroundMaker,
24 SafeMaskMaker,
25 SpectrumDatasetMaker,
26 )
27
28 # from gammapy.makers.utils import make_theta_squared_table
29 from gammapy.maps import MapAxis, RegionGeom, WcsGeom
30 from oda_api.data_products import ODAAstropyTable, PictureProduct
31 from oda_api.json import CustomJSONEncoder
32 from regions import CircleSkyRegion
33
34 hess_data = "gammapy-datasets/1.1/hess-dl3-dr1/"
35 if not (os.path.exists(hess_data)):
36 get_ipython().system("gammapy download datasets") # noqa: F821
37
38 data_store = DataStore.from_dir(hess_data)
39
40 # src_name='Crab' #http://odahub.io/ontology#AstrophysicalObject
41 # RA = 83.628700 # http://odahub.io/ontology#PointOfInterestRA
42 # DEC = 22.014700 # http://odahub.io/ontology#PointOfInterestDEC
43 src_name = "PKS 2155-304"
44 RA = 329.716938 # http://odahub.io/ontology#PointOfInterestRA
45 DEC = -30.225588 # http://odahub.io/ontology#PointOfInterestDEC
46 T1 = "2000-10-09T13:16:00.0" # http://odahub.io/ontology#StartTime
47 T2 = "2022-10-10T13:16:00.0" # http://odahub.io/ontology#EndTime
48 Radius = 2.5 # http://odahub.io/ontology#AngleDegrees
49 R_s = 0.5 # http://odahub.io/ontology#AngleDegrees
50
51 Emin = 100.0 # http://odahub.io/ontology#Energy_GeV
52 Emax = 10000.0 # http://odahub.io/ontology#Energy_GeV
53 NEbins = 20 # http://odahub.io/ontology#Integer
54
55 _galaxy_wd = os.getcwd()
56
57 with open("inputs.json", "r") as fd:
58 inp_dic = json.load(fd)
59 if "_data_product" in inp_dic.keys():
60 inp_pdic = inp_dic["_data_product"]
61 else:
62 inp_pdic = inp_dic
63
64 for vn, vv in inp_pdic.items():
65 if vn != "_selector":
66 globals()[vn] = type(globals()[vn])(vv)
67
68 T1 = Time(T1, format="isot", scale="utc").mjd
69 T2 = Time(T2, format="isot", scale="utc").mjd
70
71 dates1 = data_store.obs_table["DATE-OBS"]
72 dates2 = data_store.obs_table["DATE-END"]
73 times1 = data_store.obs_table["TIME-OBS"]
74 times2 = data_store.obs_table["TIME-END"]
75 OBSIDs = data_store.obs_table["OBS_ID"]
76 Tstart = []
77 Tstop = []
78 for i in range(len(dates1)):
79 Tstart.append(
80 Time(dates1[i] + "T" + times1[i], format="isot", scale="utc").mjd
81 )
82 Tstop.append(
83 Time(dates2[i] + "T" + times2[i], format="isot", scale="utc").mjd
84 )
85
86 RA_pnts = np.array(data_store.obs_table["RA_PNT"])
87 DEC_pnts = np.array(data_store.obs_table["DEC_PNT"])
88
89 Coords_s = SkyCoord(RA, DEC, unit="degree")
90 COORDS_pnts = SkyCoord(RA_pnts, DEC_pnts, unit="degree")
91 seps = COORDS_pnts.separation(Coords_s).deg
92
93 mask = np.where((seps < Radius) & (Tstart > T1) & (Tstop < T2))[0]
94 OBSlist = []
95 obs_ids = OBSIDs[mask]
96 if len(obs_ids) == 0:
97 message = "No data found"
98 raise RuntimeError("No data found")
99 obs_ids
100
101 observations = data_store.get_observations(obs_ids)
102
103 target_position = Coords_s
104 on_region_radius = Angle(str(R_s) + " deg")
105 on_region = CircleSkyRegion(center=target_position, radius=on_region_radius)
106 skydir = target_position.galactic
107 geom = WcsGeom.create(
108 npix=(150, 150), binsz=0.05, skydir=skydir, proj="TAN", frame="icrs"
109 )
110
111 Emin = 100.0 # http://odahub.io/ontology#Energy_GeV
112 Emax = 10000.0 # http://odahub.io/ontology#Energy_GeV
113 NEbins = 20 # http://odahub.io/ontology#Integer
114
115 energy_axis = MapAxis.from_energy_bounds(
116 Emin * 1e-3,
117 Emax * 1e-3,
118 nbin=NEbins,
119 per_decade=True,
120 unit="TeV",
121 name="energy",
122 )
123 energy_axis_true = MapAxis.from_energy_bounds(
124 0.05, 100, nbin=20, per_decade=True, unit="TeV", name="energy_true"
125 )
126
127 geom = RegionGeom.create(region=on_region, axes=[energy_axis])
128 dataset_empty = SpectrumDataset.create(
129 geom=geom, energy_axis_true=energy_axis_true
130 )
131
132 dataset_maker = SpectrumDatasetMaker(
133 containment_correction=True, selection=["counts", "exposure", "edisp"]
134 )
135 bkg_maker = ReflectedRegionsBackgroundMaker()
136 safe_mask_masker = SafeMaskMaker(methods=["aeff-max"], aeff_percent=10)
137
138 datasets = Datasets()
139
140 for obs_id, observation in zip(obs_ids, observations):
141 dataset = dataset_maker.run(
142 dataset_empty.copy(name=str(obs_id)), observation
143 )
144 dataset_on_off = bkg_maker.run(dataset, observation)
145 # dataset_on_off = safe_mask_masker.run(dataset_on_off, observation)
146 datasets.append(dataset_on_off)
147
148 print(datasets)
149
150 from pathlib import Path
151
152 path = Path("spectrum_analysis")
153 path.mkdir(exist_ok=True)
154
155 for dataset in datasets:
156 dataset.write(
157 filename=path / f"obs_{dataset.name}.fits.gz", overwrite=True
158 )
159
160 datasets = Datasets()
161
162 for obs_id in obs_ids:
163 filename = path / f"obs_{obs_id}.fits.gz"
164 datasets.append(SpectrumDatasetOnOff.read(filename))
165
166 from gammapy.modeling import Fit
167 from gammapy.modeling.models import ExpCutoffPowerLawSpectralModel, SkyModel
168
169 spectral_model = ExpCutoffPowerLawSpectralModel(
170 amplitude=1e-12 * u.Unit("cm-2 s-1 TeV-1"),
171 index=2,
172 lambda_=0.1 * u.Unit("TeV-1"),
173 reference=1 * u.TeV,
174 )
175 model = SkyModel(spectral_model=spectral_model, name="crab")
176
177 datasets.models = [model]
178
179 fit_joint = Fit()
180 result_joint = fit_joint.run(datasets=datasets)
181
182 # we make a copy here to compare it later
183 model_best_joint = model.copy()
184
185 print(result_joint)
186
187 display(result_joint.models.to_parameters_table()) # noqa: F821
188
189 e_min, e_max = Emin * 1e-3, Emax * 1e-3
190 energy_edges = np.geomspace(e_min, e_max, NEbins) * u.TeV
191
192 from gammapy.estimators import FluxPointsEstimator
193
194 fpe = FluxPointsEstimator(
195 energy_edges=energy_edges, source="crab", selection_optional="all"
196 )
197 flux_points = fpe.run(datasets=datasets)
198
199 flux_points_dataset = FluxPointsDataset(
200 data=flux_points, models=model_best_joint
201 )
202 flux_points_dataset.plot_fit()
203 # plt.show()
204 plt.savefig("Spectrum.png", format="png", bbox_inches="tight")
205
206 res = flux_points.to_table(sed_type="dnde", formatted=True)
207 np.array(res["dnde"])
208
209 bin_image = PictureProduct.from_file("Spectrum.png")
210 from astropy.table import Table
211
212 Emean = np.array(res["e_ref"])
213 Emin = np.array(res["e_min"])
214 Emax = np.array(res["e_max"])
215 flux = Emean**2 * np.array(res["dnde"])
216 flux_err = Emean**2 * np.array(res["dnde_err"])
217 data = [Emean, Emin, Emax, flux, flux_err]
218 names = (
219 "Emean[TeV]",
220 "Emin[TeV]",
221 "Emax[TeV]",
222 "Flux[TeV/cm2s]",
223 "Flux_error[TeV/cm2s]",
224 )
225 spec = ODAAstropyTable(Table(data, names=names))
226
227 picture_png = bin_image # http://odahub.io/ontology#ODAPictureProduct
228 spectrum_astropy_table = spec # http://odahub.io/ontology#ODAAstropyTable
229
230 # output gathering
231 _galaxy_meta_data = {}
232 _oda_outs = []
233 _oda_outs.append(
234 (
235 "out_Spectrum_gammapy_picture_png",
236 "picture_png_galaxy.output",
237 picture_png,
238 )
239 )
240 _oda_outs.append(
241 (
242 "out_Spectrum_gammapy_spectrum_astropy_table",
243 "spectrum_astropy_table_galaxy.output",
244 spectrum_astropy_table,
245 )
246 )
247
248 for _outn, _outfn, _outv in _oda_outs:
249 _galaxy_outfile_name = os.path.join(_galaxy_wd, _outfn)
250 if isinstance(_outv, str) and os.path.isfile(_outv):
251 shutil.move(_outv, _galaxy_outfile_name)
252 _galaxy_meta_data[_outn] = {"ext": "_sniff_"}
253 elif getattr(_outv, "write_fits_file", None):
254 _outv.write_fits_file(_galaxy_outfile_name)
255 _galaxy_meta_data[_outn] = {"ext": "fits"}
256 elif getattr(_outv, "write_file", None):
257 _outv.write_file(_galaxy_outfile_name)
258 _galaxy_meta_data[_outn] = {"ext": "_sniff_"}
259 else:
260 with open(_galaxy_outfile_name, "w") as fd:
261 json.dump(_outv, fd, cls=CustomJSONEncoder)
262 _galaxy_meta_data[_outn] = {"ext": "json"}
263
264 with open(os.path.join(_galaxy_wd, "galaxy.json"), "w") as fd:
265 json.dump(_galaxy_meta_data, fd)
266 print("*** Job finished successfully ***")