Mercurial > repos > astroteam > hess_astro_tool

comparison Spectrum.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
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)