Mercurial > repos > astroteam > hess_astro_tool

comparison Spectrum.py @ 0:02e4bb4fa10c draft

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 2991f65b25d4e6d1b69458603fce917adff40f94
author astroteam
date Mon, 19 Feb 2024 10:56:44 +0000
parents
children 593c4b45eda5
comparison
equal deleted inserted replaced
-1:000000000000 0:02e4bb4fa10c
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 matplotlib.pyplot as plt
11 import numpy as np
12 from astropy.coordinates import SkyCoord
13 from astropy.io import fits
14 from astropy.time import Time
15 from numpy import log10, sqrt
16 from oda_api.data_products import ODAAstropyTable, PictureProduct
17 from oda_api.json import CustomJSONEncoder
18
19 if os.path.exists("hess_dl3_dr1.tar.gz") == False:
20 get_ipython().system( # noqa: F821
21 "wget https://zenodo.org/record/1421099/files/hess_dl3_dr1.tar.gz"
22 )
23 get_ipython().system("tar -zxvf hess_dl3_dr1.tar.gz") # noqa: F821
24
25 # src_name='Crab' #http://odahub.io/ontology#AstrophysicalObject
26 # RA = 83.628700 # http://odahub.io/ontology#PointOfInterestRA
27 # DEC = 22.014700 # http://odahub.io/ontology#PointOfInterestDEC
28 src_name = "PKS 2155-304"
29 RA = 329.716938 # http://odahub.io/ontology#PointOfInterestRA
30 DEC = -30.225588 # http://odahub.io/ontology#PointOfInterestDEC
31
32 T1 = "2000-10-09T13:16:00.0" # http://odahub.io/ontology#StartTime
33 T2 = "2022-10-10T13:16:00.0" # http://odahub.io/ontology#EndTime
34 Radius = 2.5 # http://odahub.io/ontology#AngleDegrees
35 R_s = 0.2 # http://odahub.io/ontology#AngleDegrees
36
37 Emin = 100.0 # http://odahub.io/ontology#Energy_GeV
38 Emax = 10000.0 # http://odahub.io/ontology#Energy_GeV
39 NEbins = 20 # http://odahub.io/ontology#Integer
40
41 _galaxy_wd = os.getcwd()
42
43 with open("inputs.json", "r") as fd:
44 inp_dic = json.load(fd)
45 if "_data_product" in inp_dic.keys():
46 inp_pdic = inp_dic["_data_product"]
47 else:
48 inp_pdic = inp_dic
49
50 for vn, vv in inp_pdic.items():
51 if vn != "_selector":
52 globals()[vn] = type(globals()[vn])(vv)
53
54 Emin = Emin / 1e3
55 Emax = Emax / 1e3
56 Ebins = np.logspace(log10(Emin), log10(Emax), NEbins + 1)
57 Ebins
58 Emin = Ebins[:-1]
59 Emax = Ebins[1:]
60 Emean = sqrt(Emin * Emax)
61 lgEmean = log10(Emean)
62
63 T1 = Time(T1, format="isot", scale="utc").mjd
64 T2 = Time(T2, format="isot", scale="utc").mjd
65 message = ""
66 RA_pnts = []
67 DEC_pnts = []
68 DL3_files = []
69 OBSIDs = []
70 Tstart = []
71 Tstop = []
72 flist = os.listdir("data")
73 for f in flist:
74 if f[-7:] == "fits.gz":
75 DL3_files.append(f)
76 OBSIDs.append(int(f[20:26]))
77 hdul = fits.open("data/" + f)
78 RA_pnts.append(float(hdul[1].header["RA_PNT"]))
79 DEC_pnts.append(float(hdul[1].header["DEC_PNT"]))
80 Tstart.append(
81 Time(
82 hdul[1].header["DATE-OBS"] + "T" + hdul[1].header["TIME-OBS"],
83 format="isot",
84 scale="utc",
85 ).mjd
86 )
87 Tstop.append(
88 Time(
89 hdul[1].header["DATE-END"] + "T" + hdul[1].header["TIME-END"],
90 format="isot",
91 scale="utc",
92 ).mjd
93 )
94 hdul.close()
95
96 Coords_s = SkyCoord(RA, DEC, unit="degree")
97 COORDS_pnts = SkyCoord(RA_pnts, DEC_pnts, unit="degree")
98 seps = COORDS_pnts.separation(Coords_s).deg
99
100 mask = np.where((seps < Radius) & (Tstart > T1) & (Tstop < T2))[0]
101 OBSlist = []
102 for i in mask:
103 OBSlist.append(DL3_files[i])
104 if len(OBSlist) == 0:
105 message = "No data found"
106 raise RuntimeError("No data found")
107 message
108
109 cts_s = np.zeros(NEbins)
110 cts_b = np.zeros(NEbins)
111 Expos = np.zeros(NEbins)
112 for f in OBSlist:
113 hdul = fits.open("data/" + f)
114 RA_pnt = hdul[1].header["RA_PNT"]
115 DEC_pnt = hdul[1].header["DEC_PNT"]
116 Texp = hdul[1].header["LIVETIME"]
117 dRA = RA - RA_pnt
118 dDEC = DEC - DEC_pnt
119 RA_b = RA_pnt - dRA
120 DEC_b = DEC_pnt - dDEC
121 Coords_b = SkyCoord(RA_b, DEC_b, unit="degree")
122 Coords_pnt = SkyCoord(RA_pnt, DEC_pnt, unit="degree")
123 dist = Coords_pnt.separation(Coords_s).deg
124
125 ev = hdul["EVENTS"].data
126 ev_ra = ev["RA"]
127 ev_dec = ev["DEC"]
128 ev_en = ev["ENERGY"]
129 ev_time = ev["TIME"]
130 ev_coords = SkyCoord(ev_ra, ev_dec, unit="degree")
131 sep_s = ev_coords.separation(Coords_s).deg
132 sep_b = ev_coords.separation(Coords_b).deg
133
134 hdu = hdul["AEFF"].data
135 EEmin = hdu["ENERG_LO"][0]
136 EEmax = hdu["ENERG_HI"][0]
137 lgEE = log10(sqrt(EEmin * EEmax))
138 lgAA = log10(hdu["EFFAREA"][0] + 1e-10)
139 Thmin = hdu["THETA_LO"][0]
140 Thmax = hdu["THETA_HI"][0]
141 ind = np.argmin((Thmin - dist) ** 2)
142 Expos += 10 ** (np.interp(lgEmean, lgEE, lgAA[ind])) * Texp
143 mask = sep_s < R_s
144 cts_s += np.histogram(ev_en[mask], bins=Ebins)[0]
145 mask = sep_b < R_s
146 cts_b += np.histogram(ev_en[mask], bins=Ebins)[0]
147 hdul.close()
148
149 flux = (cts_s - cts_b) / (Emax - Emin) * Emax * Emin / (Expos * 1e4)
150 flux_err = sqrt(cts_s + cts_b) / (Emax - Emin) * Emax * Emin / (Expos * 1e4)
151 plt.errorbar(Emean, flux, yerr=flux_err, xerr=[Emean - Emin, Emax - Emean])
152 plt.xscale("log")
153 plt.yscale("log")
154 plt.xlabel("$E$, TeV")
155 plt.ylabel("$E^2 dN/dE$, erg/(cm$^2$s)")
156 plt.savefig("Spectrum.png", format="png")
157
158 bin_image = PictureProduct.from_file("Spectrum.png")
159 from astropy.table import Table
160
161 data = [Emean, Emin, Emax, flux, flux_err, cts_s, cts_b, Expos * 1e4]
162 names = (
163 "Emean[TeV]",
164 "Emin[TeV]",
165 "Emax[TeV]",
166 "Flux[TeV/cm2s]",
167 "Flux_error[TeV/cm2s]",
168 "Cts_s",
169 "Cts_b",
170 "Exposure[cm2s]",
171 )
172 spec = ODAAstropyTable(Table(data, names=names))
173
174 picture_png = bin_image # http://odahub.io/ontology#ODAPictureProduct
175 spectrum_astropy_table = spec # http://odahub.io/ontology#ODAAstropyTable
176
177 # output gathering
178 _galaxy_meta_data = {}
179 _oda_outs = []
180 _oda_outs.append(
181 ("out_Spectrum_picture_png", "picture_png_galaxy.output", picture_png)
182 )
183 _oda_outs.append(
184 (
185 "out_Spectrum_spectrum_astropy_table",
186 "spectrum_astropy_table_galaxy.output",
187 spectrum_astropy_table,
188 )
189 )
190
191 for _outn, _outfn, _outv in _oda_outs:
192 _galaxy_outfile_name = os.path.join(_galaxy_wd, _outfn)
193 if isinstance(_outv, str) and os.path.isfile(_outv):
194 shutil.move(_outv, _galaxy_outfile_name)
195 _galaxy_meta_data[_outn] = {"ext": "_sniff_"}
196 elif getattr(_outv, "write_fits_file", None):
197 _outv.write_fits_file(_galaxy_outfile_name)
198 _galaxy_meta_data[_outn] = {"ext": "fits"}
199 elif getattr(_outv, "write_file", None):
200 _outv.write_file(_galaxy_outfile_name)
201 _galaxy_meta_data[_outn] = {"ext": "_sniff_"}
202 else:
203 with open(_galaxy_outfile_name, "w") as fd:
204 json.dump(_outv, fd, cls=CustomJSONEncoder)
205 _galaxy_meta_data[_outn] = {"ext": "json"}
206
207 with open(os.path.join(_galaxy_wd, "galaxy.json"), "w") as fd:
208 json.dump(_galaxy_meta_data, fd)
209 print("*** Job finished successfully ***")