Mercurial > repos > astroteam > hess_astro_tool
changeset 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 | |
| files | Image.py Lightcurve.py Spectrum.py Spectrum_gammapy.py hess_astro_tool.xml |
| diffstat | 5 files changed, 469 insertions(+), 473 deletions(-) [+] |
line wrap: on
line diff
--- a/Image.py Mon Feb 19 10:56:44 2024 +0000 +++ b/Image.py Thu Apr 18 09:26:15 2024 +0000 @@ -14,9 +14,13 @@ from astropy.io import fits from astropy.time import Time from numpy import cos, pi +from oda_api.api import ProgressReporter from oda_api.data_products import ImageDataProduct, PictureProduct from oda_api.json import CustomJSONEncoder +pr = ProgressReporter() +pr.report_progress(stage="Progress", progress=5.0) + if os.path.exists("hess_dl3_dr1.tar.gz") == False: get_ipython().system( # noqa: F821 "wget https://zenodo.org/record/1421099/files/hess_dl3_dr1.tar.gz" @@ -28,12 +32,12 @@ DEC = 22.014700 # http://odahub.io/ontology#PointOfInterestDEC T1 = "2000-10-09T13:16:00.0" # http://odahub.io/ontology#StartTime T2 = "2022-10-10T13:16:00.0" # http://odahub.io/ontology#EndTime -Radius = 2.5 # http://odahub.io/ontology#AngleDegrees +Radius = 1.0 # http://odahub.io/ontology#AngleDegrees pixsize = ( - 0.1 # http://odahub.io/ontology#AngleDegrees ; oda:label "Pixel size" + 0.05 # http://odahub.io/ontology#AngleDegrees ; oda:label "Pixel size" ) -Emin = 100.0 # http://odahub.io/ontology#Energy_GeV -Emax = 10000.0 # http://odahub.io/ontology#Energy_GeV +Emin = 1 # http://odahub.io/ontology#Energy_TeV +Emax = 100.0 # http://odahub.io/ontology#Energy_TeV _galaxy_wd = os.getcwd() @@ -95,9 +99,14 @@ message cdec = cos(DEC * pi / 180.0) -Npix = int(4 * Radius / pixsize) + 1 -RA_bins = np.linspace(RA - Radius / cdec, RA + Radius / cdec, Npix + 1) -DEC_bins = np.linspace(DEC - Radius, DEC + Radius, Npix + 1) +Npix = int(2 * Radius / pixsize) + 1 +RA_bins = np.linspace( + RA - Npix * pixsize / cdec / 2, RA + Npix * pixsize / cdec / 2, Npix + 1 +) +DEC_bins = np.linspace( + DEC - Npix * pixsize / 2, DEC + Npix * pixsize / 2, Npix + 1 +) + image = np.zeros((Npix, Npix)) for f in OBSlist: hdul = fits.open("data/" + f) @@ -106,62 +115,76 @@ ev_dec = ev["DEC"] ev_en = ev["ENERGY"] ev_time = ev["TIME"] - h = np.histogram2d(ev_ra, ev_dec, bins=[RA_bins, DEC_bins]) + mask = (ev_en > Emin) & (ev_en < Emax) + h = np.histogram2d(ev_ra[mask], ev_dec[mask], bins=[RA_bins, DEC_bins]) image += h[0] hdul.close() +image = np.transpose(image) + plt.imshow( - np.flip(image, axis=1), - extent=(RA_bins[-1], RA_bins[0], DEC_bins[0], DEC_bins[-1]), + image, + extent=(RA_bins[0], RA_bins[-1], DEC_bins[0], DEC_bins[-1]), origin="lower", ) plt.colorbar() +plt.xlim(*plt.xlim()[::-1]) + plt.xlabel("RA, degrees") plt.ylabel("DEC,degrees") -plt.savefig("Image.png", format="png") # Create a new WCS object. The number of axes must be set # from the start w = wcs.WCS(naxis=2) -# Set up an "Airy's zenithal" projection -# Vector properties may be set with Python lists, or Numpy arrays -w.wcs.crpix = [Npix / 2.0, Npix / 2.0] -w.wcs.cdelt = np.array([pixsize / cdec, pixsize]) -w.wcs.crval = [RA, DEC] -w.wcs.ctype = ["RA---AIR", "DEC--AIR"] -w.wcs.set_pv([(2, 1, 45.0)]) +w.wcs.ctype = ["RA---CAR", "DEC--CAR"] +# we need a Plate carrée (CAR) projection since histogram is binned by ra-dec +# the peculiarity here is that CAR projection produces rectilinear grid only if CRVAL2==0 +# also, we will follow convention of RA increasing from right to left (CDELT1<0, need to flip an input image) +# otherwise, aladin-lite doesn't show it +w.wcs.crval = [RA, 0] +w.wcs.crpix = [Npix / 2.0 + 0.5, 0.5 - DEC_bins[0] / pixsize] +w.wcs.cdelt = np.array([-pixsize / cdec, pixsize]) -# Now, write out the WCS object as a FITS header header = w.to_header() -# header is an astropy.io.fits.Header object. We can use it to create a new -# PrimaryHDU and write it to a file. -hdu = fits.PrimaryHDU(image, header=header) +hdu = fits.PrimaryHDU(np.flip(image, axis=1), header=header) hdu.writeto("Image.fits", overwrite=True) hdu = fits.open("Image.fits") im = hdu[0].data -from astropy.wcs import WCS +wcs1 = wcs.WCS(hdu[0].header) +ax = plt.subplot(projection=wcs1) +plt.imshow(im, origin="lower") +plt.colorbar(label="Counts per pixel") +plt.scatter( + [RA], [DEC], marker="x", color="white", transform=ax.get_transform("world") +) +plt.text( + RA, + DEC + 0.5 * pixsize, + src_name, + color="white", + transform=ax.get_transform("world"), +) -wcs = WCS(hdu[0].header) -plt.subplot(projection=wcs) -plt.imshow(im, origin="lower") plt.grid(color="white", ls="solid") plt.xlabel("RA") plt.ylabel("Dec") +pr.report_progress(stage="Progress", progress=100.0) +plt.savefig("Image.png", format="png") bin_image = PictureProduct.from_file("Image.png") fits_image = ImageDataProduct.from_fits_file("Image.fits") -picture = bin_image # http://odahub.io/ontology#ODAPictureProduct -image = fits_image # http://odahub.io/ontology#Image +png = bin_image # http://odahub.io/ontology#ODAPictureProduct +fits = fits_image # http://odahub.io/ontology#Image # output gathering _galaxy_meta_data = {} _oda_outs = [] -_oda_outs.append(("out_Image_picture", "picture_galaxy.output", picture)) -_oda_outs.append(("out_Image_image", "image_galaxy.output", image)) +_oda_outs.append(("out_Image_png", "png_galaxy.output", png)) +_oda_outs.append(("out_Image_fits", "fits_galaxy.output", fits)) for _outn, _outfn, _outv in _oda_outs: _galaxy_outfile_name = os.path.join(_galaxy_wd, _outfn)
--- a/Lightcurve.py Mon Feb 19 10:56:44 2024 +0000 +++ b/Lightcurve.py Thu Apr 18 09:26:15 2024 +0000 @@ -21,17 +21,21 @@ "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=5.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 -T1 = "2003-10-09T13:16:00.0" # http://odahub.io/ontology#StartTime -T2 = "2005-10-10T13:16:00.0" # http://odahub.io/ontology#EndTime +T1 = "2004-11-20T13:16:00.0" # http://odahub.io/ontology#StartTime +T2 = "2004-12-20T13:16:00.0" # http://odahub.io/ontology#EndTime 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 -NTbins = 10 # http://odahub.io/ontology#Integer +Emin = 1 # http://odahub.io/ontology#Energy_TeV +Emax = 100.0 # http://odahub.io/ontology#Energy_TeV +NTbins = 30 # http://odahub.io/ontology#Integer _galaxy_wd = os.getcwd() @@ -105,7 +109,12 @@ flux_b = np.zeros(NTbins) flux_b_err = np.zeros(NTbins) Expos = np.zeros(NTbins) +src_cts = np.zeros(NTbins) +bkg_cts = np.zeros(NTbins) for count, f in enumerate(OBSlist): + pr.report_progress( + stage="Progress", progress=5.0 + 95 * count / len(OBSlist) + ) hdul = fits.open("data/" + f) RA_pnt = hdul[1].header["RA_PNT"] DEC_pnt = hdul[1].header["DEC_PNT"] @@ -124,6 +133,7 @@ ev_dec = ev["DEC"] ev_en = ev["ENERGY"] ev_time = (ev["TIME"] - Trun_start) / 86400.0 + Tbegs[count] + Nevents = len(ev) print(ev_time[0]) ev_coords = SkyCoord(ev_ra, ev_dec, unit="degree") sep_s = ev_coords.separation(Coords_s).deg @@ -138,18 +148,31 @@ Thmin = hdu["THETA_LO"][0] Thmax = hdu["THETA_HI"][0] ind = np.argmin((Thmin - dist) ** 2) - AA = AA[ind] * Texp * 1e4 - mask = np.where((sep_s < R_s)) - cts1 = np.histogram2d(ev_time[mask], ev_en[mask], bins=[Tbins, EEbins])[0] - mask = sep_b < R_s - cts2 = np.histogram2d(ev_time[mask], ev_en[mask], bins=[Tbins, EEbins])[0] - src = cts1 - cts2 - src_err = sqrt(cts1 + cts2) - flux += np.sum(src / AA, axis=1) - flux_err += np.sum(src_err / AA, axis=1) - flux_b += np.sum(cts2 / AA, axis=1) - flux_b_err += np.sum(sqrt(cts2) / AA, axis=1) + mask = EE < Emin + ind_en = len(EE[mask]) + Expos += np.histogram( + ev_time, + weights=AA[ind, ind_en] * Texp * 1e4 * np.ones(Nevents) / Nevents, + bins=Tbins, + )[0] + mask = np.where((sep_s < R_s) & (ev_en > Emin) & (ev_en < Emax)) + src_cts += np.histogram(ev_time[mask], bins=Tbins)[0] + mask = np.where((sep_b < R_s) & (ev_en > Emin) & (ev_en < Emax)) + bkg_cts += np.histogram(ev_time[mask], bins=Tbins)[0] hdul.close() +print(src_cts) +print(bkg_cts) +print(Expos) +src = src_cts - bkg_cts +src_err = sqrt(src_cts + bkg_cts) +flux = src / (Expos + 1) +flux_err = src_err / (Expos + 1) +flux_b = bkg_cts / (Expos + 1) +flux_b_err = sqrt(bkg_cts) / (Expos + 1) +# flux_err+=np.sum(src_err/AA,axis=1) +# flux_b+=np.sum(cts2/AA,axis=1) +# flux_b_err+=np.sum(sqrt(cts2)/AA,axis=1) +print(flux) if message == "": plt.errorbar( @@ -193,20 +216,14 @@ ) lc = ODAAstropyTable(Table(data, names=names)) -picture = bin_image # http://odahub.io/ontology#ODAPictureProduct -lightcurve_astropy_table = lc # http://odahub.io/ontology#ODAAstropyTable +png = bin_image # http://odahub.io/ontology#ODAPictureProduct +table = lc # http://odahub.io/ontology#ODAAstropyTable # output gathering _galaxy_meta_data = {} _oda_outs = [] -_oda_outs.append(("out_Lightcurve_picture", "picture_galaxy.output", picture)) -_oda_outs.append( - ( - "out_Lightcurve_lightcurve_astropy_table", - "lightcurve_astropy_table_galaxy.output", - lightcurve_astropy_table, - ) -) +_oda_outs.append(("out_Lightcurve_png", "png_galaxy.output", png)) +_oda_outs.append(("out_Lightcurve_table", "table_galaxy.output", table)) for _outn, _outfn, _outv in _oda_outs: _galaxy_outfile_name = os.path.join(_galaxy_wd, _outfn)
--- 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, ) )
--- a/Spectrum_gammapy.py Mon Feb 19 10:56:44 2024 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,266 +0,0 @@ -#!/usr/bin/env python -# coding: utf-8 - -# flake8: noqa - -import json -import os -import shutil - -import astropy.units as u -import matplotlib.pyplot as plt -import numpy as np -from astropy.coordinates import Angle, SkyCoord -from astropy.time import Time -from gammapy.data import DataStore -from gammapy.datasets import ( - Datasets, - FluxPointsDataset, - SpectrumDataset, - SpectrumDatasetOnOff, -) -from gammapy.makers import ( - ReflectedRegionsBackgroundMaker, - SafeMaskMaker, - SpectrumDatasetMaker, -) - -# from gammapy.makers.utils import make_theta_squared_table -from gammapy.maps import MapAxis, RegionGeom, WcsGeom -from oda_api.data_products import ODAAstropyTable, PictureProduct -from oda_api.json import CustomJSONEncoder -from regions import CircleSkyRegion - -hess_data = "gammapy-datasets/1.1/hess-dl3-dr1/" -if not (os.path.exists(hess_data)): - get_ipython().system("gammapy download datasets") # noqa: F821 - -data_store = DataStore.from_dir(hess_data) - -# 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" -RA = 329.716938 # http://odahub.io/ontology#PointOfInterestRA -DEC = -30.225588 # http://odahub.io/ontology#PointOfInterestDEC -T1 = "2000-10-09T13:16:00.0" # http://odahub.io/ontology#StartTime -T2 = "2022-10-10T13:16:00.0" # http://odahub.io/ontology#EndTime -Radius = 2.5 # http://odahub.io/ontology#AngleDegrees -R_s = 0.5 # 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 - -_galaxy_wd = os.getcwd() - -with open("inputs.json", "r") as fd: - inp_dic = json.load(fd) -if "_data_product" in inp_dic.keys(): - inp_pdic = inp_dic["_data_product"] -else: - inp_pdic = inp_dic - -for vn, vv in inp_pdic.items(): - if vn != "_selector": - globals()[vn] = type(globals()[vn])(vv) - -T1 = Time(T1, format="isot", scale="utc").mjd -T2 = Time(T2, format="isot", scale="utc").mjd - -dates1 = data_store.obs_table["DATE-OBS"] -dates2 = data_store.obs_table["DATE-END"] -times1 = data_store.obs_table["TIME-OBS"] -times2 = data_store.obs_table["TIME-END"] -OBSIDs = data_store.obs_table["OBS_ID"] -Tstart = [] -Tstop = [] -for i in range(len(dates1)): - Tstart.append( - Time(dates1[i] + "T" + times1[i], format="isot", scale="utc").mjd - ) - Tstop.append( - Time(dates2[i] + "T" + times2[i], format="isot", scale="utc").mjd - ) - -RA_pnts = np.array(data_store.obs_table["RA_PNT"]) -DEC_pnts = np.array(data_store.obs_table["DEC_PNT"]) - -Coords_s = SkyCoord(RA, DEC, unit="degree") -COORDS_pnts = SkyCoord(RA_pnts, DEC_pnts, unit="degree") -seps = COORDS_pnts.separation(Coords_s).deg - -mask = np.where((seps < Radius) & (Tstart > T1) & (Tstop < T2))[0] -OBSlist = [] -obs_ids = OBSIDs[mask] -if len(obs_ids) == 0: - message = "No data found" - raise RuntimeError("No data found") -obs_ids - -observations = data_store.get_observations(obs_ids) - -target_position = Coords_s -on_region_radius = Angle(str(R_s) + " deg") -on_region = CircleSkyRegion(center=target_position, radius=on_region_radius) -skydir = target_position.galactic -geom = WcsGeom.create( - npix=(150, 150), binsz=0.05, skydir=skydir, proj="TAN", frame="icrs" -) - -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 - -energy_axis = MapAxis.from_energy_bounds( - Emin * 1e-3, - Emax * 1e-3, - nbin=NEbins, - per_decade=True, - unit="TeV", - name="energy", -) -energy_axis_true = MapAxis.from_energy_bounds( - 0.05, 100, nbin=20, per_decade=True, unit="TeV", name="energy_true" -) - -geom = RegionGeom.create(region=on_region, axes=[energy_axis]) -dataset_empty = SpectrumDataset.create( - geom=geom, energy_axis_true=energy_axis_true -) - -dataset_maker = SpectrumDatasetMaker( - containment_correction=True, selection=["counts", "exposure", "edisp"] -) -bkg_maker = ReflectedRegionsBackgroundMaker() -safe_mask_masker = SafeMaskMaker(methods=["aeff-max"], aeff_percent=10) - -datasets = Datasets() - -for obs_id, observation in zip(obs_ids, observations): - dataset = dataset_maker.run( - dataset_empty.copy(name=str(obs_id)), observation - ) - dataset_on_off = bkg_maker.run(dataset, observation) - # dataset_on_off = safe_mask_masker.run(dataset_on_off, observation) - datasets.append(dataset_on_off) - -print(datasets) - -from pathlib import Path - -path = Path("spectrum_analysis") -path.mkdir(exist_ok=True) - -for dataset in datasets: - dataset.write( - filename=path / f"obs_{dataset.name}.fits.gz", overwrite=True - ) - -datasets = Datasets() - -for obs_id in obs_ids: - filename = path / f"obs_{obs_id}.fits.gz" - datasets.append(SpectrumDatasetOnOff.read(filename)) - -from gammapy.modeling import Fit -from gammapy.modeling.models import ExpCutoffPowerLawSpectralModel, SkyModel - -spectral_model = ExpCutoffPowerLawSpectralModel( - amplitude=1e-12 * u.Unit("cm-2 s-1 TeV-1"), - index=2, - lambda_=0.1 * u.Unit("TeV-1"), - reference=1 * u.TeV, -) -model = SkyModel(spectral_model=spectral_model, name="crab") - -datasets.models = [model] - -fit_joint = Fit() -result_joint = fit_joint.run(datasets=datasets) - -# we make a copy here to compare it later -model_best_joint = model.copy() - -print(result_joint) - -display(result_joint.models.to_parameters_table()) # noqa: F821 - -e_min, e_max = Emin * 1e-3, Emax * 1e-3 -energy_edges = np.geomspace(e_min, e_max, NEbins) * u.TeV - -from gammapy.estimators import FluxPointsEstimator - -fpe = FluxPointsEstimator( - energy_edges=energy_edges, source="crab", selection_optional="all" -) -flux_points = fpe.run(datasets=datasets) - -flux_points_dataset = FluxPointsDataset( - data=flux_points, models=model_best_joint -) -flux_points_dataset.plot_fit() -# plt.show() -plt.savefig("Spectrum.png", format="png", bbox_inches="tight") - -res = flux_points.to_table(sed_type="dnde", formatted=True) -np.array(res["dnde"]) - -bin_image = PictureProduct.from_file("Spectrum.png") -from astropy.table import Table - -Emean = np.array(res["e_ref"]) -Emin = np.array(res["e_min"]) -Emax = np.array(res["e_max"]) -flux = Emean**2 * np.array(res["dnde"]) -flux_err = Emean**2 * np.array(res["dnde_err"]) -data = [Emean, Emin, Emax, flux, flux_err] -names = ( - "Emean[TeV]", - "Emin[TeV]", - "Emax[TeV]", - "Flux[TeV/cm2s]", - "Flux_error[TeV/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 - -# output gathering -_galaxy_meta_data = {} -_oda_outs = [] -_oda_outs.append( - ( - "out_Spectrum_gammapy_picture_png", - "picture_png_galaxy.output", - picture_png, - ) -) -_oda_outs.append( - ( - "out_Spectrum_gammapy_spectrum_astropy_table", - "spectrum_astropy_table_galaxy.output", - spectrum_astropy_table, - ) -) - -for _outn, _outfn, _outv in _oda_outs: - _galaxy_outfile_name = os.path.join(_galaxy_wd, _outfn) - if isinstance(_outv, str) and os.path.isfile(_outv): - shutil.move(_outv, _galaxy_outfile_name) - _galaxy_meta_data[_outn] = {"ext": "_sniff_"} - elif getattr(_outv, "write_fits_file", None): - _outv.write_fits_file(_galaxy_outfile_name) - _galaxy_meta_data[_outn] = {"ext": "fits"} - elif getattr(_outv, "write_file", None): - _outv.write_file(_galaxy_outfile_name) - _galaxy_meta_data[_outn] = {"ext": "_sniff_"} - else: - with open(_galaxy_outfile_name, "w") as fd: - json.dump(_outv, fd, cls=CustomJSONEncoder) - _galaxy_meta_data[_outn] = {"ext": "json"} - -with open(os.path.join(_galaxy_wd, "galaxy.json"), "w") as fd: - json.dump(_galaxy_meta_data, fd) -print("*** Job finished successfully ***")
--- a/hess_astro_tool.xml Mon Feb 19 10:56:44 2024 +0000 +++ b/hess_astro_tool.xml Thu Apr 18 09:26:15 2024 +0000 @@ -1,102 +1,88 @@ -<tool id="hess_astro_tool" name="HESS" version="0.0.1+galaxy0" profile="23.0"> +<tool id="hess_astro_tool" name="HESS" version="0.0.2+galaxy0" profile="23.0"> <requirements> - <requirement type="package" version="8.21.0">ipython</requirement> + <requirement type="package" version="8.22.2">ipython</requirement> <requirement type="package" version="5.3.4">astropy</requirement> <requirement type="package" version="1.11.4">scipy</requirement> - <requirement type="package" version="1.10.13">pydantic</requirement> + <requirement type="package" version="1.10.14">pydantic</requirement> <requirement type="package" version="4.66.2">tqdm</requirement> <requirement type="package" version="1.1">gammapy</requirement> - <requirement type="package" version="3.8.3">matplotlib</requirement> - <requirement type="package" version="1.2.12">oda-api</requirement> + <requirement type="package" version="3.8.4">matplotlib</requirement> + <requirement type="package" version="1.2.15">oda-api</requirement> <!--Requirements string 'nb2workflow[cwl,service,rdf,mmoda]>=1.3.30 ' can't be converted automatically. Please add the galaxy/conda requirement manually or modify the requirements file!--> - <requirement type="package" version="7.16.0">nbconvert</requirement> - <requirement type="package" version="1.20.3">wget</requirement> + <requirement type="package" version="7.16.3">nbconvert</requirement> + <requirement type="package" version="1.21.4">wget</requirement> </requirements> <command detect_errors="exit_code">ipython '$__tool_directory__/${_data_product._selector}.py'</command> <configfiles> - <inputs name="inputs" filename="inputs.json" /> + <inputs name="inputs" filename="inputs.json" data_style="paths" /> </configfiles> <inputs> <conditional name="_data_product"> <param name="_selector" type="select" label="Data Product"> <option value="Image" selected="true">Image</option> <option value="Spectrum" selected="false">Spectrum</option> - <option value="Spectrum_gammapy" selected="false">Spectrum_gammapy</option> <option value="Lightcurve" selected="false">Lightcurve</option> </param> <when value="Image"> <param name="src_name" type="text" value="Crab" label="src_name" /> - <param name="RA" type="float" value="83.6287" label="RA Units: deg" /> - <param name="DEC" type="float" value="22.0147" label="DEC Units: deg" /> + <param name="RA" type="float" value="83.6287" label="RA (unit: deg)" /> + <param name="DEC" type="float" value="22.0147" label="DEC (unit: deg)" /> <param name="T1" type="text" value="2000-10-09T13:16:00.0" label="T1" /> <param name="T2" type="text" value="2022-10-10T13:16:00.0" label="T2" /> - <param name="Radius" type="float" value="2.5" label="Radius Units: deg" /> - <param name="pixsize" type="float" value="0.1" label="Pixel size Units: deg" /> - <param name="Emin" type="float" value="100.0" label="Emin Units: GeV" /> - <param name="Emax" type="float" value="10000.0" label="Emax Units: GeV" /> + <param name="Radius" type="float" value="1.0" label="Radius (unit: deg)" /> + <param name="pixsize" type="float" value="0.05" label="Pixel size (unit: deg)" /> + <param name="Emin" type="integer" value="1" label="Emin (unit: TeV)" /> + <param name="Emax" type="float" value="100.0" label="Emax (unit: TeV)" /> </when> <when value="Spectrum"> <param name="src_name" type="text" value="PKS 2155-304" label="src_name" /> - <param name="RA" type="float" value="329.716938" label="RA Units: deg" /> - <param name="DEC" type="float" value="-30.225588" label="DEC Units: deg" /> + <param name="RA" type="float" value="329.716938" label="RA (unit: deg)" /> + <param name="DEC" type="float" value="-30.225588" label="DEC (unit: deg)" /> <param name="T1" type="text" value="2000-10-09T13:16:00.0" label="T1" /> <param name="T2" type="text" value="2022-10-10T13:16:00.0" label="T2" /> - <param name="Radius" type="float" value="2.5" label="Radius Units: deg" /> - <param name="R_s" type="float" value="0.2" label="R_s Units: deg" /> - <param name="Emin" type="float" value="100.0" label="Emin Units: GeV" /> - <param name="Emax" type="float" value="10000.0" label="Emax Units: GeV" /> - <param name="NEbins" type="integer" value="20" label="NEbins" /> - </when> - <when value="Spectrum_gammapy"> - <param name="src_name" type="text" value="PKS 2155-304" label="src_name" /> - <param name="RA" type="float" value="329.716938" label="RA Units: deg" /> - <param name="DEC" type="float" value="-30.225588" label="DEC Units: deg" /> - <param name="T1" type="text" value="2000-10-09T13:16:00.0" label="T1" /> - <param name="T2" type="text" value="2022-10-10T13:16:00.0" label="T2" /> - <param name="Radius" type="float" value="2.5" label="Radius Units: deg" /> - <param name="R_s" type="float" value="0.5" label="R_s Units: deg" /> - <param name="Emin" type="float" value="100.0" label="Emin Units: GeV" /> - <param name="Emax" type="float" value="10000.0" label="Emax Units: GeV" /> - <param name="NEbins" type="integer" value="20" label="NEbins" /> + <param name="Radius" type="float" value="2.5" label="Radius (unit: deg)" /> + <param name="R_s" type="float" value="0.2" label="R_s (unit: deg)" /> + <param name="Emin" type="float" value="0.1" label="Emin (unit: TeV)" /> + <param name="Emax" type="float" value="100.0" label="Emax (unit: TeV)" /> + <param name="NEbins" type="integer" value="30" label="NEbins" /> + <param name="Efit_min" type="float" value="0.2" label="Efit_min (unit: TeV)" /> + <param name="Efit_max" type="float" value="10.0" label="Efit_max (unit: TeV)" /> </when> <when value="Lightcurve"> <param name="src_name" type="text" value="Crab" label="src_name" /> - <param name="RA" type="float" value="83.6287" label="RA Units: deg" /> - <param name="DEC" type="float" value="22.0147" label="DEC Units: deg" /> - <param name="T1" type="text" value="2003-10-09T13:16:00.0" label="T1" /> - <param name="T2" type="text" value="2005-10-10T13:16:00.0" label="T2" /> - <param name="Radius" type="float" value="2.5" label="Radius Units: deg" /> - <param name="R_s" type="float" value="0.2" label="R_s Units: deg" /> - <param name="Emin" type="float" value="100.0" label="Emin Units: GeV" /> - <param name="Emax" type="float" value="10000.0" label="Emax Units: GeV" /> - <param name="NTbins" type="integer" value="10" label="NTbins" /> + <param name="RA" type="float" value="83.6287" label="RA (unit: deg)" /> + <param name="DEC" type="float" value="22.0147" label="DEC (unit: deg)" /> + <param name="T1" type="text" value="2004-11-20T13:16:00.0" label="T1" /> + <param name="T2" type="text" value="2004-12-20T13:16:00.0" label="T2" /> + <param name="Radius" type="float" value="2.5" label="Radius (unit: deg)" /> + <param name="R_s" type="float" value="0.2" label="R_s (unit: deg)" /> + <param name="Emin" type="integer" value="1" label="Emin (unit: TeV)" /> + <param name="Emax" type="float" value="100.0" label="Emax (unit: TeV)" /> + <param name="NTbins" type="integer" value="30" label="NTbins" /> </when> </conditional> </inputs> <outputs> - <data label="${tool.name} -> Image picture" name="out_Image_picture" format="auto" from_work_dir="picture_galaxy.output"> + <data label="${tool.name} -> Image png" name="out_Image_png" format="auto" from_work_dir="png_galaxy.output"> <filter>_data_product['_selector'] == 'Image'</filter> </data> - <data label="${tool.name} -> Image image" name="out_Image_image" format="auto" from_work_dir="image_galaxy.output"> + <data label="${tool.name} -> Image fits" name="out_Image_fits" format="auto" from_work_dir="fits_galaxy.output"> <filter>_data_product['_selector'] == 'Image'</filter> </data> - <data label="${tool.name} -> Spectrum picture_png" name="out_Spectrum_picture_png" format="auto" from_work_dir="picture_png_galaxy.output"> + <data label="${tool.name} -> Spectrum png" name="out_Spectrum_png" format="auto" from_work_dir="png_galaxy.output"> <filter>_data_product['_selector'] == 'Spectrum'</filter> </data> - <data label="${tool.name} -> Spectrum spectrum_astropy_table" name="out_Spectrum_spectrum_astropy_table" format="auto" from_work_dir="spectrum_astropy_table_galaxy.output"> + <data label="${tool.name} -> Spectrum table_confidence_contour" name="out_Spectrum_table_confidence_contour" format="auto" from_work_dir="table_confidence_contour_galaxy.output"> <filter>_data_product['_selector'] == 'Spectrum'</filter> </data> - <data label="${tool.name} -> Spectrum_gammapy picture_png" name="out_Spectrum_gammapy_picture_png" format="auto" from_work_dir="picture_png_galaxy.output"> - <filter>_data_product['_selector'] == 'Spectrum_gammapy'</filter> + <data label="${tool.name} -> Spectrum table_spectrum" name="out_Spectrum_table_spectrum" format="auto" from_work_dir="table_spectrum_galaxy.output"> + <filter>_data_product['_selector'] == 'Spectrum'</filter> </data> - <data label="${tool.name} -> Spectrum_gammapy spectrum_astropy_table" name="out_Spectrum_gammapy_spectrum_astropy_table" format="auto" from_work_dir="spectrum_astropy_table_galaxy.output"> - <filter>_data_product['_selector'] == 'Spectrum_gammapy'</filter> - </data> - <data label="${tool.name} -> Lightcurve picture" name="out_Lightcurve_picture" format="auto" from_work_dir="picture_galaxy.output"> + <data label="${tool.name} -> Lightcurve png" name="out_Lightcurve_png" format="auto" from_work_dir="png_galaxy.output"> <filter>_data_product['_selector'] == 'Lightcurve'</filter> </data> - <data label="${tool.name} -> Lightcurve lightcurve_astropy_table" name="out_Lightcurve_lightcurve_astropy_table" format="auto" from_work_dir="lightcurve_astropy_table_galaxy.output"> + <data label="${tool.name} -> Lightcurve table" name="out_Lightcurve_table" format="auto" from_work_dir="table_galaxy.output"> <filter>_data_product['_selector'] == 'Lightcurve'</filter> </data> </outputs> @@ -109,16 +95,16 @@ <param name="DEC" value="22.0147" /> <param name="T1" value="2000-10-09T13:16:00.0" /> <param name="T2" value="2022-10-10T13:16:00.0" /> - <param name="Radius" value="2.5" /> - <param name="pixsize" value="0.1" /> - <param name="Emin" value="100.0" /> - <param name="Emax" value="10000.0" /> + <param name="Radius" value="1.0" /> + <param name="pixsize" value="0.05" /> + <param name="Emin" value="1" /> + <param name="Emax" value="100.0" /> </conditional> <assert_stdout> <has_text text="*** Job finished successfully ***" /> </assert_stdout> </test> - <test expect_num_outputs="2"> + <test expect_num_outputs="3"> <conditional name="_data_product"> <param name="_selector" value="Spectrum" /> <param name="src_name" value="PKS 2155-304" /> @@ -128,27 +114,11 @@ <param name="T2" value="2022-10-10T13:16:00.0" /> <param name="Radius" value="2.5" /> <param name="R_s" value="0.2" /> - <param name="Emin" value="100.0" /> - <param name="Emax" value="10000.0" /> - <param name="NEbins" value="20" /> - </conditional> - <assert_stdout> - <has_text text="*** Job finished successfully ***" /> - </assert_stdout> - </test> - <test expect_num_outputs="2"> - <conditional name="_data_product"> - <param name="_selector" value="Spectrum_gammapy" /> - <param name="src_name" value="PKS 2155-304" /> - <param name="RA" value="329.716938" /> - <param name="DEC" value="-30.225588" /> - <param name="T1" value="2000-10-09T13:16:00.0" /> - <param name="T2" value="2022-10-10T13:16:00.0" /> - <param name="Radius" value="2.5" /> - <param name="R_s" value="0.5" /> - <param name="Emin" value="100.0" /> - <param name="Emax" value="10000.0" /> - <param name="NEbins" value="20" /> + <param name="Emin" value="0.1" /> + <param name="Emax" value="100.0" /> + <param name="NEbins" value="30" /> + <param name="Efit_min" value="0.2" /> + <param name="Efit_max" value="10.0" /> </conditional> <assert_stdout> <has_text text="*** Job finished successfully ***" /> @@ -160,13 +130,13 @@ <param name="src_name" value="Crab" /> <param name="RA" value="83.6287" /> <param name="DEC" value="22.0147" /> - <param name="T1" value="2003-10-09T13:16:00.0" /> - <param name="T2" value="2005-10-10T13:16:00.0" /> + <param name="T1" value="2004-11-20T13:16:00.0" /> + <param name="T2" value="2004-12-20T13:16:00.0" /> <param name="Radius" value="2.5" /> <param name="R_s" value="0.2" /> - <param name="Emin" value="100.0" /> - <param name="Emax" value="10000.0" /> - <param name="NTbins" value="10" /> + <param name="Emin" value="1" /> + <param name="Emax" value="100.0" /> + <param name="NTbins" value="30" /> </conditional> <assert_stdout> <has_text text="*** Job finished successfully ***" /> @@ -185,31 +155,37 @@ the time range that can be adjusted setting the ``T1`` (start time) and ``T2`` (stop time) parameters. +For the spectra, the analysis performs histogramming of events in +energy, in the number ``NEbins`` of energy bins homogeneously spaces in +logarithm of energy, beterrn ``Emin`` and ``Emax``. Conversion of the +counts to the physical flux units is done by dividing by the exposure +time and effective area that is extracted from the Instrument Response +Functions (IRF). The result is shown with black data points. This simple +estimate of flux in energy bins does not take into account the event +energy estimation errors. + +The source signal is extracted from a circular region of the radius +``R_s``\ (adjustable parameter) around the source posiiton. The +background estimate is done using the “wobble” method, from a position +opposite with respect to the signal extraction region in the camera, +with respect to the center of the telescope field-of-view. + +A powerlaw fit to the spectrum is done using forward folding method, +properly taking into account the error of energy estimation. The results +of the spectral fitting include the 68% confidence contour for the +spectral parameters (slope and normalisaiton of the powerlaw spectrum). +The first entry in the “contour” table is the best-fit value of the fit. +The energy range of fitting can be adjusted with the ``Efit_min``, +``Efit_max`` parameters. + The lightcurves are produced by hystrogramming of the events in time, in the number ``NTbins`` of time intervals of equalt time width between ``T1`` (start time) and ``T2`` (stop time). The events are selected in a -desired energy range between ``Emin`` and ``Emax`` from a circular -region of the radius ``R_s`` (in degrees) around the source position -``RA``,\ ``DEC``. Conversion of the counts to the physical flux units is -done by dividing by the exposure time and effective area that is -extracted from the Instrument Response Functions (IRF). - -For the spectra, two alternative tools are considered. The service -``Spectrum`` performs histogramming of the events in energy, in the -number ``NEbins`` of energy bins homogeneously spaces in logarithm of -energy, beterrn ``Emin`` and ``Emax``. Conversion of the counts to the -physical flux units is done by dividing by the exposure time and -effective area that is extracted from the IRF. This method does not take -into account the energy bias and can result in a wrong spectral shape at -the low-energy threshold where the bias is strongest. - -An alternative spectral extraction is done using -`Gammapy <https://gammapy.org/>`__, following the script `Spectrum -Analysis <https://docs.gammapy.org/0.18.2/tutorials/spectrum_analysis.html>`__. -It performs forward folding of the spectral model (a cut-off powerlaw by -default) with the IRF and fits the folded model to the binned count data -in the energy range between ``Emin`` and ``Emax`` in ``NEbins`` -logarithmically spaced energy bins. +desired energy range between ``Emin`` and ``Emax``. The source and +backgorund counts are computed in the same way as for the spectral +fitting. Conversion of the counts to the physical flux units is done by +dividing by the exposure time and effective area that is extracted from +the IRF. </help> <citations> <citation type="doi">10.5281/zenodo.1421098</citation>