Mercurial > repos > astroteam > hess_astro_tool
diff 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 |
line wrap: on
line diff
--- a/Spectrum.py Mon Feb 19 10:56:44 2024 +0000 +++ b/Spectrum.py Thu Apr 18 09:26:15 2024 +0000 @@ -21,11 +21,15 @@ "wget https://zenodo.org/record/1421099/files/hess_dl3_dr1.tar.gz" ) get_ipython().system("tar -zxvf hess_dl3_dr1.tar.gz") # noqa: F821 +from oda_api.api import ProgressReporter + +pr = ProgressReporter() +pr.report_progress(stage="Progress", progress=0.0) # src_name='Crab' #http://odahub.io/ontology#AstrophysicalObject # RA = 83.628700 # http://odahub.io/ontology#PointOfInterestRA # DEC = 22.014700 # http://odahub.io/ontology#PointOfInterestDEC -src_name = "PKS 2155-304" +src_name = "PKS 2155-304" # http://odahub.io/ontology#AstrophysicalObject RA = 329.716938 # http://odahub.io/ontology#PointOfInterestRA DEC = -30.225588 # http://odahub.io/ontology#PointOfInterestDEC @@ -34,9 +38,12 @@ Radius = 2.5 # http://odahub.io/ontology#AngleDegrees R_s = 0.2 # http://odahub.io/ontology#AngleDegrees -Emin = 100.0 # http://odahub.io/ontology#Energy_GeV -Emax = 10000.0 # http://odahub.io/ontology#Energy_GeV -NEbins = 20 # http://odahub.io/ontology#Integer +Emin = 0.1 # http://odahub.io/ontology#Energy_TeV +Emax = 100.0 # http://odahub.io/ontology#Energy_TeV +NEbins = 30 # http://odahub.io/ontology#Integer + +Efit_min = 0.2 # http://odahub.io/ontology#Energy_TeV +Efit_max = 10.0 # http://odahub.io/ontology#Energy_TeV _galaxy_wd = os.getcwd() @@ -51,14 +58,22 @@ if vn != "_selector": globals()[vn] = type(globals()[vn])(vv) -Emin = Emin / 1e3 -Emax = Emax / 1e3 +E0 = 1.0 + +def model_dNdE(E, N, Gam): + return N * (E / E0) ** (Gam) + +def model_rate(E1, E2, N, Gam): + dEE = E2 - E1 + EE = sqrt(E1 * E2) + return model_dNdE(EE, N, Gam) * dEE + Ebins = np.logspace(log10(Emin), log10(Emax), NEbins + 1) -Ebins -Emin = Ebins[:-1] -Emax = Ebins[1:] -Emean = sqrt(Emin * Emax) -lgEmean = log10(Emean) +Emins = Ebins[:-1] +Emaxs = Ebins[1:] +Emeans = sqrt(Emins * Emaxs) +lgEmeans = log10(Emeans) +dE = Ebins[1:] - Ebins[:-1] T1 = Time(T1, format="isot", scale="utc").mjd T2 = Time(T2, format="isot", scale="utc").mjd @@ -104,16 +119,35 @@ if len(OBSlist) == 0: message = "No data found" raise RuntimeError("No data found") -message +offaxis = seps[mask] +Tstart = np.array(Tstart)[mask] +print("Found", len(Tstart), "pointings") -cts_s = np.zeros(NEbins) -cts_b = np.zeros(NEbins) -Expos = np.zeros(NEbins) -for f in OBSlist: - hdul = fits.open("data/" + f) +ind = 0 +pointing = OBSlist[ind] +hdul = fits.open("data/" + pointing) +RMF = hdul["EDISP"].data +ENERG_LO = RMF["ENERG_LO"][0] +ENERG_HI = RMF["ENERG_HI"][0] +MIGRA_LO = RMF["MIGRA_LO"][0] # MIGRA_bins=np.linspace(0.2,5,161) +MIGRA_HI = RMF["MIGRA_HI"][0] +MIGRA = (MIGRA_LO + MIGRA_HI) / 2.0 +ENERG = sqrt(ENERG_LO * ENERG_HI) +dENERG = ENERG_HI - ENERG_LO + +cts_s = [] +cts_b = [] +Eff_area = [] +Eff_area_interp = [] +Texp = [] +RMFs = [] +for ind in range(len(OBSlist)): + pointing = OBSlist[ind] + hdul = fits.open("data/" + pointing) + RA_pnt = hdul[1].header["RA_PNT"] DEC_pnt = hdul[1].header["DEC_PNT"] - Texp = hdul[1].header["LIVETIME"] + Texp.append(hdul[1].header["LIVETIME"]) dRA = RA - RA_pnt dDEC = DEC - DEC_pnt RA_b = RA_pnt - dRA @@ -122,6 +156,26 @@ Coords_pnt = SkyCoord(RA_pnt, DEC_pnt, unit="degree") dist = Coords_pnt.separation(Coords_s).deg + RMF = hdul["EDISP"].data + mask = RMF["THETA_LO"] < dist + ind_th = len(RMF["THETA_LO"][mask]) - 1 + RMF_th = RMF["MATRIX"][0][ind_th] + RMF_interp = np.zeros((len(Emeans), len(ENERG))) + for k in range(len(ENERG)): + dp_dErec = RMF_th[:, k] / ENERG[k] + Erec = MIGRA * ENERG[k] + dp_dErec_interp = ( + np.interp(Emeans, Erec, dp_dErec) + * (Emeans > min(Erec)) + * (Emeans < max(Erec)) + ) + RMF_interp[:, k] = dp_dErec_interp + RMFs.append(RMF_interp) + + AEFF = hdul["AEFF"].data + Eff_area.append(AEFF["EFFAREA"][0][ind_th]) + Eff_area_interp.append(np.interp(Emeans, ENERG, Eff_area[-1])) + ev = hdul["EVENTS"].data ev_ra = ev["RA"] ev_dec = ev["DEC"] @@ -131,60 +185,252 @@ sep_s = ev_coords.separation(Coords_s).deg sep_b = ev_coords.separation(Coords_b).deg - hdu = hdul["AEFF"].data - EEmin = hdu["ENERG_LO"][0] - EEmax = hdu["ENERG_HI"][0] - lgEE = log10(sqrt(EEmin * EEmax)) - lgAA = log10(hdu["EFFAREA"][0] + 1e-10) - Thmin = hdu["THETA_LO"][0] - Thmax = hdu["THETA_HI"][0] - ind = np.argmin((Thmin - dist) ** 2) - Expos += 10 ** (np.interp(lgEmean, lgEE, lgAA[ind])) * Texp mask = sep_s < R_s - cts_s += np.histogram(ev_en[mask], bins=Ebins)[0] + cts_s.append(np.histogram(ev_en[mask], bins=Ebins)[0]) mask = sep_b < R_s - cts_b += np.histogram(ev_en[mask], bins=Ebins)[0] + cts_b.append(np.histogram(ev_en[mask], bins=Ebins)[0]) hdul.close() -flux = (cts_s - cts_b) / (Emax - Emin) * Emax * Emin / (Expos * 1e4) -flux_err = sqrt(cts_s + cts_b) / (Emax - Emin) * Emax * Emin / (Expos * 1e4) -plt.errorbar(Emean, flux, yerr=flux_err, xerr=[Emean - Emin, Emax - Emean]) +cts_s = np.array(cts_s) +cts_b = np.array(cts_b) +Eff_area = np.array(Eff_area) +Eff_area_interp = np.array(Eff_area_interp) +Texp = np.array(Texp) + +cts_s_tot = sum(cts_s) +cts_b_tot = sum(cts_b) +src_tot = cts_s_tot - cts_b_tot +src_tot_err = sqrt(cts_s_tot + cts_b_tot) +Expos_tot_interp = sum(Eff_area_interp * np.outer(Texp, np.ones(NEbins))) * 1e4 +flux_tot = src_tot / (Expos_tot_interp + 1) / (Emaxs - Emins) * Emaxs * Emins +flux_tot_err = ( + src_tot_err / (Expos_tot_interp + 1) / (Emaxs - Emins) * Emaxs * Emins +) +print( + "Total source counts:", + sum(cts_s_tot), + "; background counts", + sum(cts_b_tot), +) + +def model_cts_Erec(N, Gam): + model_ENERG = model_rate(ENERG_LO, ENERG_HI, N, Gam) + res = np.zeros(NEbins) + for ind in range(len(OBSlist)): + model_counts_ENERG = model_ENERG * Eff_area[ind] * 1e4 * Texp[ind] + for k in range(len(ENERG)): + res += model_counts_ENERG[k] * RMFs[ind][:, k] * dE + return res + +def chi2(p): + N, slope = p + counts = model_cts_Erec(N, slope) + m = Emeans > Efit_min + m &= Emeans < Efit_max + m &= src_tot_err > 0.0 + chi2 = (((counts[m] - src_tot[m]) / src_tot_err[m]) ** 2).sum() + dof = sum(m) + # print(N,slope,chi2) + return chi2, dof + +plt.errorbar( + Emeans, + src_tot, + src_tot_err, + xerr=[Emeans - Emins, Emaxs - Emeans], + linestyle="none", +) + +plt.axvline(Efit_min, color="black") +plt.axvline(Efit_max, color="black") plt.xscale("log") plt.yscale("log") +N = 4e-11 +Gam = -2.7 +plt.plot(Emeans, model_cts_Erec(N, Gam)) +chi2([N, Gam])[0] + +# 1) find 90% confidence contour scanning over a wide parameter space +Norm_max = 1e-10 +Norm_min = 1e-12 +Norm_bins = 100 +Gam_min = -1.0 +Gam_max = -5.0 +Gam_bins = 100 +Ns = np.linspace(Norm_min, Norm_max, Norm_bins) +Gams = np.linspace(Gam_min, Gam_max, Gam_bins) +chi2_map = np.zeros((Norm_bins, Gam_bins)) +Norm_best = Norm_min +Gam_best = Gam_min +chi2_best = 1e10 +for i, N in enumerate(Ns): + pr.report_progress(stage="Progress", progress=5 + 45 * i / Norm_bins) + for j, Gam in enumerate(Gams): + chi2_map[i, j] = chi2([N, Gam])[0] + if chi2_map[i, j] < chi2_best: + Norm_best = N + Gam_best = Gam + chi2_best = chi2_map[i, j] +print(Norm_best, Gam_best) +# 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])) + +# 68% contour from https://ui.adsabs.harvard.edu/abs/1976ApJ...208..177L/abstract for two-parameter fit +# 90% contour from https://ui.adsabs.harvard.edu/abs/1976ApJ...208..177L/abstract for two-parameter fit + +cnt = plt.contour( + Gams, Ns, chi2_map, levels=[np.amin(chi2_map) + 4.61], colors="red" +) +plt.scatter([Gam_best], [Norm_best], marker="x", color="red") +# plt.colorbar() +print(np.amin(chi2_map)) + +cont = cnt.get_paths()[0].vertices +gammas = cont[:, 0] +norms = cont[:, 1] + +# 2) refine with 68% contour calculation within the initial 90% countour +Norm_max = max(1.1 * norms) +Norm_min = min(0.9 * norms) +Norm_bins = 50 +Gam_min = min(gammas - 0.2) +Gam_max = max(gammas + 0.2) +Gam_bins = 50 +Ns = np.linspace(Norm_min, Norm_max, Norm_bins) +Gams = np.linspace(Gam_min, Gam_max, Gam_bins) +chi2_map = np.zeros((Norm_bins, Gam_bins)) +Norm_best = Norm_min +Gam_best = Gam_min +chi2_best = 1e10 +for i, N in enumerate(Ns): + pr.report_progress(stage="Progress", progress=50 + 50 * i / Norm_bins) + for j, Gam in enumerate(Gams): + chi2_map[i, j] = chi2([N, Gam])[0] + if chi2_map[i, j] < chi2_best: + Norm_best = N + Gam_best = Gam + chi2_best = chi2_map[i, j] +print(Norm_best, Gam_best) +# 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])) + +# 68% contour from https://ui.adsabs.harvard.edu/abs/1976ApJ...208..177L/abstract for two-parameter fit +cnt = plt.contour( + Gams, Ns, chi2_map, levels=[np.amin(chi2_map) + 2.3], colors="black" +) +cont = cnt.get_paths()[0].vertices +gammas = cont[:, 0] +norms = cont[:, 1] + +plt.scatter([Gam_best], [Norm_best], marker="x", color="black") +# plt.colorbar() +print(chi2([Norm_best, Gam_best])) +plt.xlabel("Slope") +plt.ylabel(r"Norm, 1/(TeV cm$^2$ s)") +plt.savefig("Contour.png", format="png", bbox_inches="tight") + +plt.figure(figsize=(8, 6)) +x = np.logspace(log10(Efit_min), log10(Efit_max), 10) +ymax = np.zeros(10) +ymin = np.ones(10) +for i in range(len(gammas)): + y = model_dNdE(x, norms[i], gammas[i]) * x**2 + ymax = np.maximum(y, ymax) + ymin = np.minimum(y, ymin) + # plt.plot(x,y) +plt.fill_between(x, ymin, ymax, alpha=0.2) +plt.plot(x, model_dNdE(x, Norm_best, Gam_best) * x**2) + +plt.errorbar( + Emeans, + flux_tot, + flux_tot_err, + xerr=[Emeans - Emins, Emaxs - Emeans], + linestyle="none", + color="black", + linewidth=2, +) + +# plt.axvspan(0, Efit_min, alpha=0.2, color='black') +# plt.axvspan(Efit_max, 1000, alpha=0.2, color='black') + +plt.xscale("log") +plt.yscale("log") +plt.xlim(0.1, 100) plt.xlabel("$E$, TeV") -plt.ylabel("$E^2 dN/dE$, erg/(cm$^2$s)") -plt.savefig("Spectrum.png", format="png") +plt.ylabel("$E^2 dN/dE$, TeV/(cm$^2$ s)") +plt.savefig("Spectrum.png", format="png", bbox_inches="tight") + +print("Here") +new_hdul = fits.HDUList() +form = str(len(ENERG)) + "E" +for i in range(len(OBSlist)): + col1 = fits.Column(name="COUNTS", format="E", array=cts_s[i]) + col2 = fits.Column(name="BACKGROUND", format="E", array=cts_b[i]) + hdu = fits.BinTableHDU.from_columns([col1, col2], name="COUNTS") + new_hdul.append(hdu) + col = fits.Column(format=form, array=RMFs[i], name="RMF") + hdu = fits.BinTableHDU.from_columns([col], name="RMF") + new_hdul.append(hdu) + col = fits.Column(format="E", array=Eff_area[i], name="ARF") + hdu = fits.BinTableHDU.from_columns([col], name="ARF") + new_hdul.append(hdu) + +col1 = fits.Column(name="E_MIN", format="E", unit="TeV", array=Emins) +col2 = fits.Column(name="E_MAX", format="E", unit="TeV", array=Emaxs) +cols = fits.ColDefs([col1, col2]) +hdu = fits.BinTableHDU.from_columns(cols, name="Erec_BOUNDS") +new_hdul.append(hdu) +col1 = fits.Column(name="ENERG_LO", format="E", unit="TeV", array=ENERG_LO) +col2 = fits.Column(name="ENERG_HI", format="E", unit="TeV", array=ENERG_HI) +cols = fits.ColDefs([col1, col2]) +hdu = fits.BinTableHDU.from_columns(cols, name="Etrue_BOUNDS") +new_hdul.append(hdu) +# new_hdul.writeto('Spectra.fits',overwrite=True) + +print("and here") bin_image = PictureProduct.from_file("Spectrum.png") +bin_image1 = PictureProduct.from_file("Contour.png") from astropy.table import Table -data = [Emean, Emin, Emax, flux, flux_err, cts_s, cts_b, Expos * 1e4] +data = [Emeans, Emins, Emaxs, flux_tot, flux_tot_err] names = ( "Emean[TeV]", "Emin[TeV]", "Emax[TeV]", "Flux[TeV/cm2s]", "Flux_error[TeV/cm2s]", - "Cts_s", - "Cts_b", - "Exposure[cm2s]", ) spec = ODAAstropyTable(Table(data, names=names)) -picture_png = bin_image # http://odahub.io/ontology#ODAPictureProduct -spectrum_astropy_table = spec # http://odahub.io/ontology#ODAAstropyTable +data = [ + np.concatenate(([Gam_best], gammas)), + np.concatenate(([Norm_best], norms)), +] +names = ["Gamma", "Norm_1TeV[1/(TeV cm2 s)]"] +error_ellipse = ODAAstropyTable(Table(data, names=names)) + +png = bin_image # http://odahub.io/ontology#ODAPictureProduct +table_confidence_contour = ( + error_ellipse # http://odahub.io/ontology#ODAAstropyTable +) +table_spectrum = spec # http://odahub.io/ontology#ODAAstropyTable # output gathering _galaxy_meta_data = {} _oda_outs = [] +_oda_outs.append(("out_Spectrum_png", "png_galaxy.output", png)) _oda_outs.append( - ("out_Spectrum_picture_png", "picture_png_galaxy.output", picture_png) + ( + "out_Spectrum_table_confidence_contour", + "table_confidence_contour_galaxy.output", + table_confidence_contour, + ) ) _oda_outs.append( ( - "out_Spectrum_spectrum_astropy_table", - "spectrum_astropy_table_galaxy.output", - spectrum_astropy_table, + "out_Spectrum_table_spectrum", + "table_spectrum_galaxy.output", + table_spectrum, ) )