Mercurial > repos > astroteam > hess_astro_tool
comparison Spectrum.py @ 1:593c4b45eda5 draft default tip
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 |
---|---|
19 if os.path.exists("hess_dl3_dr1.tar.gz") == False: | 19 if os.path.exists("hess_dl3_dr1.tar.gz") == False: |
20 get_ipython().system( # noqa: F821 | 20 get_ipython().system( # noqa: F821 |
21 "wget https://zenodo.org/record/1421099/files/hess_dl3_dr1.tar.gz" | 21 "wget https://zenodo.org/record/1421099/files/hess_dl3_dr1.tar.gz" |
22 ) | 22 ) |
23 get_ipython().system("tar -zxvf hess_dl3_dr1.tar.gz") # noqa: F821 | 23 get_ipython().system("tar -zxvf hess_dl3_dr1.tar.gz") # noqa: F821 |
24 from oda_api.api import ProgressReporter | |
25 | |
26 pr = ProgressReporter() | |
27 pr.report_progress(stage="Progress", progress=0.0) | |
24 | 28 |
25 # src_name='Crab' #http://odahub.io/ontology#AstrophysicalObject | 29 # src_name='Crab' #http://odahub.io/ontology#AstrophysicalObject |
26 # RA = 83.628700 # http://odahub.io/ontology#PointOfInterestRA | 30 # RA = 83.628700 # http://odahub.io/ontology#PointOfInterestRA |
27 # DEC = 22.014700 # http://odahub.io/ontology#PointOfInterestDEC | 31 # DEC = 22.014700 # http://odahub.io/ontology#PointOfInterestDEC |
28 src_name = "PKS 2155-304" | 32 src_name = "PKS 2155-304" # http://odahub.io/ontology#AstrophysicalObject |
29 RA = 329.716938 # http://odahub.io/ontology#PointOfInterestRA | 33 RA = 329.716938 # http://odahub.io/ontology#PointOfInterestRA |
30 DEC = -30.225588 # http://odahub.io/ontology#PointOfInterestDEC | 34 DEC = -30.225588 # http://odahub.io/ontology#PointOfInterestDEC |
31 | 35 |
32 T1 = "2000-10-09T13:16:00.0" # http://odahub.io/ontology#StartTime | 36 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 | 37 T2 = "2022-10-10T13:16:00.0" # http://odahub.io/ontology#EndTime |
34 Radius = 2.5 # http://odahub.io/ontology#AngleDegrees | 38 Radius = 2.5 # http://odahub.io/ontology#AngleDegrees |
35 R_s = 0.2 # http://odahub.io/ontology#AngleDegrees | 39 R_s = 0.2 # http://odahub.io/ontology#AngleDegrees |
36 | 40 |
37 Emin = 100.0 # http://odahub.io/ontology#Energy_GeV | 41 Emin = 0.1 # http://odahub.io/ontology#Energy_TeV |
38 Emax = 10000.0 # http://odahub.io/ontology#Energy_GeV | 42 Emax = 100.0 # http://odahub.io/ontology#Energy_TeV |
39 NEbins = 20 # http://odahub.io/ontology#Integer | 43 NEbins = 30 # http://odahub.io/ontology#Integer |
44 | |
45 Efit_min = 0.2 # http://odahub.io/ontology#Energy_TeV | |
46 Efit_max = 10.0 # http://odahub.io/ontology#Energy_TeV | |
40 | 47 |
41 _galaxy_wd = os.getcwd() | 48 _galaxy_wd = os.getcwd() |
42 | 49 |
43 with open("inputs.json", "r") as fd: | 50 with open("inputs.json", "r") as fd: |
44 inp_dic = json.load(fd) | 51 inp_dic = json.load(fd) |
49 | 56 |
50 for vn, vv in inp_pdic.items(): | 57 for vn, vv in inp_pdic.items(): |
51 if vn != "_selector": | 58 if vn != "_selector": |
52 globals()[vn] = type(globals()[vn])(vv) | 59 globals()[vn] = type(globals()[vn])(vv) |
53 | 60 |
54 Emin = Emin / 1e3 | 61 E0 = 1.0 |
55 Emax = Emax / 1e3 | 62 |
63 def model_dNdE(E, N, Gam): | |
64 return N * (E / E0) ** (Gam) | |
65 | |
66 def model_rate(E1, E2, N, Gam): | |
67 dEE = E2 - E1 | |
68 EE = sqrt(E1 * E2) | |
69 return model_dNdE(EE, N, Gam) * dEE | |
70 | |
56 Ebins = np.logspace(log10(Emin), log10(Emax), NEbins + 1) | 71 Ebins = np.logspace(log10(Emin), log10(Emax), NEbins + 1) |
57 Ebins | 72 Emins = Ebins[:-1] |
58 Emin = Ebins[:-1] | 73 Emaxs = Ebins[1:] |
59 Emax = Ebins[1:] | 74 Emeans = sqrt(Emins * Emaxs) |
60 Emean = sqrt(Emin * Emax) | 75 lgEmeans = log10(Emeans) |
61 lgEmean = log10(Emean) | 76 dE = Ebins[1:] - Ebins[:-1] |
62 | 77 |
63 T1 = Time(T1, format="isot", scale="utc").mjd | 78 T1 = Time(T1, format="isot", scale="utc").mjd |
64 T2 = Time(T2, format="isot", scale="utc").mjd | 79 T2 = Time(T2, format="isot", scale="utc").mjd |
65 message = "" | 80 message = "" |
66 RA_pnts = [] | 81 RA_pnts = [] |
102 for i in mask: | 117 for i in mask: |
103 OBSlist.append(DL3_files[i]) | 118 OBSlist.append(DL3_files[i]) |
104 if len(OBSlist) == 0: | 119 if len(OBSlist) == 0: |
105 message = "No data found" | 120 message = "No data found" |
106 raise RuntimeError("No data found") | 121 raise RuntimeError("No data found") |
107 message | 122 offaxis = seps[mask] |
108 | 123 Tstart = np.array(Tstart)[mask] |
109 cts_s = np.zeros(NEbins) | 124 print("Found", len(Tstart), "pointings") |
110 cts_b = np.zeros(NEbins) | 125 |
111 Expos = np.zeros(NEbins) | 126 ind = 0 |
112 for f in OBSlist: | 127 pointing = OBSlist[ind] |
113 hdul = fits.open("data/" + f) | 128 hdul = fits.open("data/" + pointing) |
129 RMF = hdul["EDISP"].data | |
130 ENERG_LO = RMF["ENERG_LO"][0] | |
131 ENERG_HI = RMF["ENERG_HI"][0] | |
132 MIGRA_LO = RMF["MIGRA_LO"][0] # MIGRA_bins=np.linspace(0.2,5,161) | |
133 MIGRA_HI = RMF["MIGRA_HI"][0] | |
134 MIGRA = (MIGRA_LO + MIGRA_HI) / 2.0 | |
135 ENERG = sqrt(ENERG_LO * ENERG_HI) | |
136 dENERG = ENERG_HI - ENERG_LO | |
137 | |
138 cts_s = [] | |
139 cts_b = [] | |
140 Eff_area = [] | |
141 Eff_area_interp = [] | |
142 Texp = [] | |
143 RMFs = [] | |
144 for ind in range(len(OBSlist)): | |
145 pointing = OBSlist[ind] | |
146 hdul = fits.open("data/" + pointing) | |
147 | |
114 RA_pnt = hdul[1].header["RA_PNT"] | 148 RA_pnt = hdul[1].header["RA_PNT"] |
115 DEC_pnt = hdul[1].header["DEC_PNT"] | 149 DEC_pnt = hdul[1].header["DEC_PNT"] |
116 Texp = hdul[1].header["LIVETIME"] | 150 Texp.append(hdul[1].header["LIVETIME"]) |
117 dRA = RA - RA_pnt | 151 dRA = RA - RA_pnt |
118 dDEC = DEC - DEC_pnt | 152 dDEC = DEC - DEC_pnt |
119 RA_b = RA_pnt - dRA | 153 RA_b = RA_pnt - dRA |
120 DEC_b = DEC_pnt - dDEC | 154 DEC_b = DEC_pnt - dDEC |
121 Coords_b = SkyCoord(RA_b, DEC_b, unit="degree") | 155 Coords_b = SkyCoord(RA_b, DEC_b, unit="degree") |
122 Coords_pnt = SkyCoord(RA_pnt, DEC_pnt, unit="degree") | 156 Coords_pnt = SkyCoord(RA_pnt, DEC_pnt, unit="degree") |
123 dist = Coords_pnt.separation(Coords_s).deg | 157 dist = Coords_pnt.separation(Coords_s).deg |
158 | |
159 RMF = hdul["EDISP"].data | |
160 mask = RMF["THETA_LO"] < dist | |
161 ind_th = len(RMF["THETA_LO"][mask]) - 1 | |
162 RMF_th = RMF["MATRIX"][0][ind_th] | |
163 RMF_interp = np.zeros((len(Emeans), len(ENERG))) | |
164 for k in range(len(ENERG)): | |
165 dp_dErec = RMF_th[:, k] / ENERG[k] | |
166 Erec = MIGRA * ENERG[k] | |
167 dp_dErec_interp = ( | |
168 np.interp(Emeans, Erec, dp_dErec) | |
169 * (Emeans > min(Erec)) | |
170 * (Emeans < max(Erec)) | |
171 ) | |
172 RMF_interp[:, k] = dp_dErec_interp | |
173 RMFs.append(RMF_interp) | |
174 | |
175 AEFF = hdul["AEFF"].data | |
176 Eff_area.append(AEFF["EFFAREA"][0][ind_th]) | |
177 Eff_area_interp.append(np.interp(Emeans, ENERG, Eff_area[-1])) | |
124 | 178 |
125 ev = hdul["EVENTS"].data | 179 ev = hdul["EVENTS"].data |
126 ev_ra = ev["RA"] | 180 ev_ra = ev["RA"] |
127 ev_dec = ev["DEC"] | 181 ev_dec = ev["DEC"] |
128 ev_en = ev["ENERGY"] | 182 ev_en = ev["ENERGY"] |
129 ev_time = ev["TIME"] | 183 ev_time = ev["TIME"] |
130 ev_coords = SkyCoord(ev_ra, ev_dec, unit="degree") | 184 ev_coords = SkyCoord(ev_ra, ev_dec, unit="degree") |
131 sep_s = ev_coords.separation(Coords_s).deg | 185 sep_s = ev_coords.separation(Coords_s).deg |
132 sep_b = ev_coords.separation(Coords_b).deg | 186 sep_b = ev_coords.separation(Coords_b).deg |
133 | 187 |
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 | 188 mask = sep_s < R_s |
144 cts_s += np.histogram(ev_en[mask], bins=Ebins)[0] | 189 cts_s.append(np.histogram(ev_en[mask], bins=Ebins)[0]) |
145 mask = sep_b < R_s | 190 mask = sep_b < R_s |
146 cts_b += np.histogram(ev_en[mask], bins=Ebins)[0] | 191 cts_b.append(np.histogram(ev_en[mask], bins=Ebins)[0]) |
147 hdul.close() | 192 hdul.close() |
148 | 193 |
149 flux = (cts_s - cts_b) / (Emax - Emin) * Emax * Emin / (Expos * 1e4) | 194 cts_s = np.array(cts_s) |
150 flux_err = sqrt(cts_s + cts_b) / (Emax - Emin) * Emax * Emin / (Expos * 1e4) | 195 cts_b = np.array(cts_b) |
151 plt.errorbar(Emean, flux, yerr=flux_err, xerr=[Emean - Emin, Emax - Emean]) | 196 Eff_area = np.array(Eff_area) |
197 Eff_area_interp = np.array(Eff_area_interp) | |
198 Texp = np.array(Texp) | |
199 | |
200 cts_s_tot = sum(cts_s) | |
201 cts_b_tot = sum(cts_b) | |
202 src_tot = cts_s_tot - cts_b_tot | |
203 src_tot_err = sqrt(cts_s_tot + cts_b_tot) | |
204 Expos_tot_interp = sum(Eff_area_interp * np.outer(Texp, np.ones(NEbins))) * 1e4 | |
205 flux_tot = src_tot / (Expos_tot_interp + 1) / (Emaxs - Emins) * Emaxs * Emins | |
206 flux_tot_err = ( | |
207 src_tot_err / (Expos_tot_interp + 1) / (Emaxs - Emins) * Emaxs * Emins | |
208 ) | |
209 print( | |
210 "Total source counts:", | |
211 sum(cts_s_tot), | |
212 "; background counts", | |
213 sum(cts_b_tot), | |
214 ) | |
215 | |
216 def model_cts_Erec(N, Gam): | |
217 model_ENERG = model_rate(ENERG_LO, ENERG_HI, N, Gam) | |
218 res = np.zeros(NEbins) | |
219 for ind in range(len(OBSlist)): | |
220 model_counts_ENERG = model_ENERG * Eff_area[ind] * 1e4 * Texp[ind] | |
221 for k in range(len(ENERG)): | |
222 res += model_counts_ENERG[k] * RMFs[ind][:, k] * dE | |
223 return res | |
224 | |
225 def chi2(p): | |
226 N, slope = p | |
227 counts = model_cts_Erec(N, slope) | |
228 m = Emeans > Efit_min | |
229 m &= Emeans < Efit_max | |
230 m &= src_tot_err > 0.0 | |
231 chi2 = (((counts[m] - src_tot[m]) / src_tot_err[m]) ** 2).sum() | |
232 dof = sum(m) | |
233 # print(N,slope,chi2) | |
234 return chi2, dof | |
235 | |
236 plt.errorbar( | |
237 Emeans, | |
238 src_tot, | |
239 src_tot_err, | |
240 xerr=[Emeans - Emins, Emaxs - Emeans], | |
241 linestyle="none", | |
242 ) | |
243 | |
244 plt.axvline(Efit_min, color="black") | |
245 plt.axvline(Efit_max, color="black") | |
152 plt.xscale("log") | 246 plt.xscale("log") |
153 plt.yscale("log") | 247 plt.yscale("log") |
248 N = 4e-11 | |
249 Gam = -2.7 | |
250 plt.plot(Emeans, model_cts_Erec(N, Gam)) | |
251 chi2([N, Gam])[0] | |
252 | |
253 # 1) find 90% confidence contour scanning over a wide parameter space | |
254 Norm_max = 1e-10 | |
255 Norm_min = 1e-12 | |
256 Norm_bins = 100 | |
257 Gam_min = -1.0 | |
258 Gam_max = -5.0 | |
259 Gam_bins = 100 | |
260 Ns = np.linspace(Norm_min, Norm_max, Norm_bins) | |
261 Gams = np.linspace(Gam_min, Gam_max, Gam_bins) | |
262 chi2_map = np.zeros((Norm_bins, Gam_bins)) | |
263 Norm_best = Norm_min | |
264 Gam_best = Gam_min | |
265 chi2_best = 1e10 | |
266 for i, N in enumerate(Ns): | |
267 pr.report_progress(stage="Progress", progress=5 + 45 * i / Norm_bins) | |
268 for j, Gam in enumerate(Gams): | |
269 chi2_map[i, j] = chi2([N, Gam])[0] | |
270 if chi2_map[i, j] < chi2_best: | |
271 Norm_best = N | |
272 Gam_best = Gam | |
273 chi2_best = chi2_map[i, j] | |
274 print(Norm_best, Gam_best) | |
275 # plt.imshow(chi2_map,vmax=np.amin(chi2_map)+4,origin='lower',extent=[Gams[0],Gams[-1],Ns[0],Ns[-1]],aspect=(Gams[-1]-Gams[0])/(Ns[-1]-Ns[0])) | |
276 | |
277 # 68% contour from https://ui.adsabs.harvard.edu/abs/1976ApJ...208..177L/abstract for two-parameter fit | |
278 # 90% contour from https://ui.adsabs.harvard.edu/abs/1976ApJ...208..177L/abstract for two-parameter fit | |
279 | |
280 cnt = plt.contour( | |
281 Gams, Ns, chi2_map, levels=[np.amin(chi2_map) + 4.61], colors="red" | |
282 ) | |
283 plt.scatter([Gam_best], [Norm_best], marker="x", color="red") | |
284 # plt.colorbar() | |
285 print(np.amin(chi2_map)) | |
286 | |
287 cont = cnt.get_paths()[0].vertices | |
288 gammas = cont[:, 0] | |
289 norms = cont[:, 1] | |
290 | |
291 # 2) refine with 68% contour calculation within the initial 90% countour | |
292 Norm_max = max(1.1 * norms) | |
293 Norm_min = min(0.9 * norms) | |
294 Norm_bins = 50 | |
295 Gam_min = min(gammas - 0.2) | |
296 Gam_max = max(gammas + 0.2) | |
297 Gam_bins = 50 | |
298 Ns = np.linspace(Norm_min, Norm_max, Norm_bins) | |
299 Gams = np.linspace(Gam_min, Gam_max, Gam_bins) | |
300 chi2_map = np.zeros((Norm_bins, Gam_bins)) | |
301 Norm_best = Norm_min | |
302 Gam_best = Gam_min | |
303 chi2_best = 1e10 | |
304 for i, N in enumerate(Ns): | |
305 pr.report_progress(stage="Progress", progress=50 + 50 * i / Norm_bins) | |
306 for j, Gam in enumerate(Gams): | |
307 chi2_map[i, j] = chi2([N, Gam])[0] | |
308 if chi2_map[i, j] < chi2_best: | |
309 Norm_best = N | |
310 Gam_best = Gam | |
311 chi2_best = chi2_map[i, j] | |
312 print(Norm_best, Gam_best) | |
313 # plt.imshow(chi2_map,vmax=np.amin(chi2_map)+4,origin='lower',extent=[Gams[0],Gams[-1],Ns[0],Ns[-1]],aspect=(Gams[-1]-Gams[0])/(Ns[-1]-Ns[0])) | |
314 | |
315 # 68% contour from https://ui.adsabs.harvard.edu/abs/1976ApJ...208..177L/abstract for two-parameter fit | |
316 cnt = plt.contour( | |
317 Gams, Ns, chi2_map, levels=[np.amin(chi2_map) + 2.3], colors="black" | |
318 ) | |
319 cont = cnt.get_paths()[0].vertices | |
320 gammas = cont[:, 0] | |
321 norms = cont[:, 1] | |
322 | |
323 plt.scatter([Gam_best], [Norm_best], marker="x", color="black") | |
324 # plt.colorbar() | |
325 print(chi2([Norm_best, Gam_best])) | |
326 plt.xlabel("Slope") | |
327 plt.ylabel(r"Norm, 1/(TeV cm$^2$ s)") | |
328 plt.savefig("Contour.png", format="png", bbox_inches="tight") | |
329 | |
330 plt.figure(figsize=(8, 6)) | |
331 x = np.logspace(log10(Efit_min), log10(Efit_max), 10) | |
332 ymax = np.zeros(10) | |
333 ymin = np.ones(10) | |
334 for i in range(len(gammas)): | |
335 y = model_dNdE(x, norms[i], gammas[i]) * x**2 | |
336 ymax = np.maximum(y, ymax) | |
337 ymin = np.minimum(y, ymin) | |
338 # plt.plot(x,y) | |
339 plt.fill_between(x, ymin, ymax, alpha=0.2) | |
340 plt.plot(x, model_dNdE(x, Norm_best, Gam_best) * x**2) | |
341 | |
342 plt.errorbar( | |
343 Emeans, | |
344 flux_tot, | |
345 flux_tot_err, | |
346 xerr=[Emeans - Emins, Emaxs - Emeans], | |
347 linestyle="none", | |
348 color="black", | |
349 linewidth=2, | |
350 ) | |
351 | |
352 # plt.axvspan(0, Efit_min, alpha=0.2, color='black') | |
353 # plt.axvspan(Efit_max, 1000, alpha=0.2, color='black') | |
354 | |
355 plt.xscale("log") | |
356 plt.yscale("log") | |
357 plt.xlim(0.1, 100) | |
154 plt.xlabel("$E$, TeV") | 358 plt.xlabel("$E$, TeV") |
155 plt.ylabel("$E^2 dN/dE$, erg/(cm$^2$s)") | 359 plt.ylabel("$E^2 dN/dE$, TeV/(cm$^2$ s)") |
156 plt.savefig("Spectrum.png", format="png") | 360 plt.savefig("Spectrum.png", format="png", bbox_inches="tight") |
157 | 361 |
362 print("Here") | |
363 new_hdul = fits.HDUList() | |
364 form = str(len(ENERG)) + "E" | |
365 | |
366 for i in range(len(OBSlist)): | |
367 col1 = fits.Column(name="COUNTS", format="E", array=cts_s[i]) | |
368 col2 = fits.Column(name="BACKGROUND", format="E", array=cts_b[i]) | |
369 hdu = fits.BinTableHDU.from_columns([col1, col2], name="COUNTS") | |
370 new_hdul.append(hdu) | |
371 col = fits.Column(format=form, array=RMFs[i], name="RMF") | |
372 hdu = fits.BinTableHDU.from_columns([col], name="RMF") | |
373 new_hdul.append(hdu) | |
374 col = fits.Column(format="E", array=Eff_area[i], name="ARF") | |
375 hdu = fits.BinTableHDU.from_columns([col], name="ARF") | |
376 new_hdul.append(hdu) | |
377 | |
378 col1 = fits.Column(name="E_MIN", format="E", unit="TeV", array=Emins) | |
379 col2 = fits.Column(name="E_MAX", format="E", unit="TeV", array=Emaxs) | |
380 cols = fits.ColDefs([col1, col2]) | |
381 hdu = fits.BinTableHDU.from_columns(cols, name="Erec_BOUNDS") | |
382 new_hdul.append(hdu) | |
383 col1 = fits.Column(name="ENERG_LO", format="E", unit="TeV", array=ENERG_LO) | |
384 col2 = fits.Column(name="ENERG_HI", format="E", unit="TeV", array=ENERG_HI) | |
385 cols = fits.ColDefs([col1, col2]) | |
386 hdu = fits.BinTableHDU.from_columns(cols, name="Etrue_BOUNDS") | |
387 new_hdul.append(hdu) | |
388 # new_hdul.writeto('Spectra.fits',overwrite=True) | |
389 | |
390 print("and here") | |
158 bin_image = PictureProduct.from_file("Spectrum.png") | 391 bin_image = PictureProduct.from_file("Spectrum.png") |
392 bin_image1 = PictureProduct.from_file("Contour.png") | |
159 from astropy.table import Table | 393 from astropy.table import Table |
160 | 394 |
161 data = [Emean, Emin, Emax, flux, flux_err, cts_s, cts_b, Expos * 1e4] | 395 data = [Emeans, Emins, Emaxs, flux_tot, flux_tot_err] |
162 names = ( | 396 names = ( |
163 "Emean[TeV]", | 397 "Emean[TeV]", |
164 "Emin[TeV]", | 398 "Emin[TeV]", |
165 "Emax[TeV]", | 399 "Emax[TeV]", |
166 "Flux[TeV/cm2s]", | 400 "Flux[TeV/cm2s]", |
167 "Flux_error[TeV/cm2s]", | 401 "Flux_error[TeV/cm2s]", |
168 "Cts_s", | |
169 "Cts_b", | |
170 "Exposure[cm2s]", | |
171 ) | 402 ) |
172 spec = ODAAstropyTable(Table(data, names=names)) | 403 spec = ODAAstropyTable(Table(data, names=names)) |
173 | 404 |
174 picture_png = bin_image # http://odahub.io/ontology#ODAPictureProduct | 405 data = [ |
175 spectrum_astropy_table = spec # http://odahub.io/ontology#ODAAstropyTable | 406 np.concatenate(([Gam_best], gammas)), |
407 np.concatenate(([Norm_best], norms)), | |
408 ] | |
409 names = ["Gamma", "Norm_1TeV[1/(TeV cm2 s)]"] | |
410 error_ellipse = ODAAstropyTable(Table(data, names=names)) | |
411 | |
412 png = bin_image # http://odahub.io/ontology#ODAPictureProduct | |
413 table_confidence_contour = ( | |
414 error_ellipse # http://odahub.io/ontology#ODAAstropyTable | |
415 ) | |
416 table_spectrum = spec # http://odahub.io/ontology#ODAAstropyTable | |
176 | 417 |
177 # output gathering | 418 # output gathering |
178 _galaxy_meta_data = {} | 419 _galaxy_meta_data = {} |
179 _oda_outs = [] | 420 _oda_outs = [] |
180 _oda_outs.append( | 421 _oda_outs.append(("out_Spectrum_png", "png_galaxy.output", png)) |
181 ("out_Spectrum_picture_png", "picture_png_galaxy.output", picture_png) | |
182 ) | |
183 _oda_outs.append( | 422 _oda_outs.append( |
184 ( | 423 ( |
185 "out_Spectrum_spectrum_astropy_table", | 424 "out_Spectrum_table_confidence_contour", |
186 "spectrum_astropy_table_galaxy.output", | 425 "table_confidence_contour_galaxy.output", |
187 spectrum_astropy_table, | 426 table_confidence_contour, |
427 ) | |
428 ) | |
429 _oda_outs.append( | |
430 ( | |
431 "out_Spectrum_table_spectrum", | |
432 "table_spectrum_galaxy.output", | |
433 table_spectrum, | |
188 ) | 434 ) |
189 ) | 435 ) |
190 | 436 |
191 for _outn, _outfn, _outv in _oda_outs: | 437 for _outn, _outfn, _outv in _oda_outs: |
192 _galaxy_outfile_name = os.path.join(_galaxy_wd, _outfn) | 438 _galaxy_outfile_name = os.path.join(_galaxy_wd, _outfn) |