Mercurial > repos > astroteam > hess_astro_tool
comparison Lightcurve.py @ 0:02e4bb4fa10c draft
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 ***") |