Mercurial > repos > astroteam > hess_astro_tool

comparison Lightcurve.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 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 T1 = "2003-10-09T13:16:00.0" # http://odahub.io/ontology#StartTime
29 T2 = "2005-10-10T13:16:00.0" # http://odahub.io/ontology#EndTime
30 Radius = 2.5 # http://odahub.io/ontology#AngleDegrees
31 R_s = 0.2 # http://odahub.io/ontology#AngleDegrees
32 Emin = 100.0 # http://odahub.io/ontology#Energy_GeV
33 Emax = 10000.0 # http://odahub.io/ontology#Energy_GeV
34 NTbins = 10 # http://odahub.io/ontology#Integer
35
36 _galaxy_wd = os.getcwd()
37
38 with open("inputs.json", "r") as fd:
39 inp_dic = json.load(fd)
40 if "_data_product" in inp_dic.keys():
41 inp_pdic = inp_dic["_data_product"]
42 else:
43 inp_pdic = inp_dic
44
45 for vn, vv in inp_pdic.items():
46 if vn != "_selector":
47 globals()[vn] = type(globals()[vn])(vv)
48
49 T1 = Time(T1, format="isot", scale="utc").mjd
50 T2 = Time(T2, format="isot", scale="utc").mjd
51 message = ""
52 RA_pnts = []
53 DEC_pnts = []
54 DL3_files = []
55 OBSIDs = []
56 Tstart = []
57 Tstop = []
58 flist = os.listdir("data")
59 for f in flist:
60 if f[-7:] == "fits.gz":
61 DL3_files.append(f)
62 OBSIDs.append(int(f[20:26]))
63 hdul = fits.open("data/" + f)
64 RA_pnts.append(float(hdul[1].header["RA_PNT"]))
65 DEC_pnts.append(float(hdul[1].header["DEC_PNT"]))
66 Tstart.append(
67 Time(
68 hdul[1].header["DATE-OBS"] + "T" + hdul[1].header["TIME-OBS"],
69 format="isot",
70 scale="utc",
71 ).mjd
72 )
73 Tstop.append(
74 Time(
75 hdul[1].header["DATE-END"] + "T" + hdul[1].header["TIME-END"],
76 format="isot",
77 scale="utc",
78 ).mjd
79 )
80 hdul.close()
81
82 Coords_s = SkyCoord(RA, DEC, unit="degree")
83 COORDS_pnts = SkyCoord(RA_pnts, DEC_pnts, unit="degree")
84 seps = COORDS_pnts.separation(Coords_s).deg
85
86 mask = np.where((seps < Radius) & (Tstart > T1) & (Tstop < T2))[0]
87 OBSlist = []
88 Tbegs = []
89 for i in mask:
90 OBSlist.append(DL3_files[i])
91 Tbegs.append(Tstart[i])
92 if len(OBSlist) == 0:
93 message = "No data found"
94 raise RuntimeError("No data found")
95 message
96
97 Tbins = np.linspace(T1, T2, NTbins + 1)
98 Tmin = Tbins[:-1]
99 Tmax = Tbins[1:]
100 Tmean = (Tmin + Tmax) / 2.0
101 Tbins
102
103 flux = np.zeros(NTbins)
104 flux_err = np.zeros(NTbins)
105 flux_b = np.zeros(NTbins)
106 flux_b_err = np.zeros(NTbins)
107 Expos = np.zeros(NTbins)
108 for count, f in enumerate(OBSlist):
109 hdul = fits.open("data/" + f)
110 RA_pnt = hdul[1].header["RA_PNT"]
111 DEC_pnt = hdul[1].header["DEC_PNT"]
112 Texp = hdul[1].header["LIVETIME"]
113 Trun_start = hdul[1].header["TSTART"]
114 dRA = RA - RA_pnt
115 dDEC = DEC - DEC_pnt
116 RA_b = RA_pnt - dRA
117 DEC_b = DEC_pnt - dDEC
118 Coords_b = SkyCoord(RA_b, DEC_b, unit="degree")
119 Coords_pnt = SkyCoord(RA_pnt, DEC_pnt, unit="degree")
120 dist = Coords_pnt.separation(Coords_s).deg
121
122 ev = hdul["EVENTS"].data
123 ev_ra = ev["RA"]
124 ev_dec = ev["DEC"]
125 ev_en = ev["ENERGY"]
126 ev_time = (ev["TIME"] - Trun_start) / 86400.0 + Tbegs[count]
127 print(ev_time[0])
128 ev_coords = SkyCoord(ev_ra, ev_dec, unit="degree")
129 sep_s = ev_coords.separation(Coords_s).deg
130 sep_b = ev_coords.separation(Coords_b).deg
131
132 hdu = hdul["AEFF"].data
133 EEmin = hdu["ENERG_LO"][0]
134 EEmax = hdu["ENERG_HI"][0]
135 EE = sqrt(EEmin * EEmax)
136 EEbins = np.concatenate((EEmin, [EEmax[-1]]))
137 AA = hdu["EFFAREA"][0] + 1e-10
138 Thmin = hdu["THETA_LO"][0]
139 Thmax = hdu["THETA_HI"][0]
140 ind = np.argmin((Thmin - dist) ** 2)
141 AA = AA[ind] * Texp * 1e4
142 mask = np.where((sep_s < R_s))
143 cts1 = np.histogram2d(ev_time[mask], ev_en[mask], bins=[Tbins, EEbins])[0]
144 mask = sep_b < R_s
145 cts2 = np.histogram2d(ev_time[mask], ev_en[mask], bins=[Tbins, EEbins])[0]
146 src = cts1 - cts2
147 src_err = sqrt(cts1 + cts2)
148 flux += np.sum(src / AA, axis=1)
149 flux_err += np.sum(src_err / AA, axis=1)
150 flux_b += np.sum(cts2 / AA, axis=1)
151 flux_b_err += np.sum(sqrt(cts2) / AA, axis=1)
152 hdul.close()
153
154 if message == "":
155 plt.errorbar(
156 Tmean,
157 flux,
158 yerr=flux_err,
159 xerr=[Tmean - Tmin, Tmax - Tmean],
160 linestyle="none",
161 label="source",
162 )
163 plt.errorbar(
164 Tmean,
165 flux_b,
166 yerr=flux_b_err,
167 xerr=[Tmean - Tmin, Tmax - Tmean],
168 linestyle="none",
169 label="background",
170 )
171 plt.xlabel("Time, MJD")
172 plt.ylabel("Flux, cts/cm$^2$s")
173 plt.yscale("log")
174 ymin = min(min(flux - flux_err), min(flux_b - flux_b_err))
175 ymax = max(max(flux + flux_err), max(flux_b + flux_b_err))
176 plt.ylim(ymin / 2.0, 2 * ymax)
177 plt.legend(loc="lower left")
178 plt.savefig("Lightcurve.png", format="png")
179
180 if message == "":
181 bin_image = PictureProduct.from_file("Lightcurve.png")
182 from astropy.table import Table
183
184 data = [Tmean, Tmin, Tmax, flux, flux_err, flux_b, flux_b_err]
185 names = (
186 "Tmean[MJD]",
187 "Tmin[MJD]",
188 "Tmax[MJD]",
189 "Flux[counts/cm2s]",
190 "Flux_error[counts/cm2s]",
191 "Background[counts/cm2s]",
192 "Background_error[counts/cm2s]",
193 )
194 lc = ODAAstropyTable(Table(data, names=names))
195
196 picture = bin_image # http://odahub.io/ontology#ODAPictureProduct
197 lightcurve_astropy_table = lc # http://odahub.io/ontology#ODAAstropyTable
198
199 # output gathering
200 _galaxy_meta_data = {}
201 _oda_outs = []
202 _oda_outs.append(("out_Lightcurve_picture", "picture_galaxy.output", picture))
203 _oda_outs.append(
204 (
205 "out_Lightcurve_lightcurve_astropy_table",
206 "lightcurve_astropy_table_galaxy.output",
207 lightcurve_astropy_table,
208 )
209 )
210
211 for _outn, _outfn, _outv in _oda_outs:
212 _galaxy_outfile_name = os.path.join(_galaxy_wd, _outfn)
213 if isinstance(_outv, str) and os.path.isfile(_outv):
214 shutil.move(_outv, _galaxy_outfile_name)
215 _galaxy_meta_data[_outn] = {"ext": "_sniff_"}
216 elif getattr(_outv, "write_fits_file", None):
217 _outv.write_fits_file(_galaxy_outfile_name)
218 _galaxy_meta_data[_outn] = {"ext": "fits"}
219 elif getattr(_outv, "write_file", None):
220 _outv.write_file(_galaxy_outfile_name)
221 _galaxy_meta_data[_outn] = {"ext": "_sniff_"}
222 else:
223 with open(_galaxy_outfile_name, "w") as fd:
224 json.dump(_outv, fd, cls=CustomJSONEncoder)
225 _galaxy_meta_data[_outn] = {"ext": "json"}
226
227 with open(os.path.join(_galaxy_wd, "galaxy.json"), "w") as fd:
228 json.dump(_galaxy_meta_data, fd)
229 print("*** Job finished successfully ***")