Mercurial > repos > astroteam > hess_astro_tool
diff Lightcurve.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/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)