Mercurial > repos > astroteam > hess_astro_tool
comparison 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 |
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=5.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 T1 = "2003-10-09T13:16:00.0" # http://odahub.io/ontology#StartTime | 32 T1 = "2004-11-20T13:16:00.0" # http://odahub.io/ontology#StartTime |
29 T2 = "2005-10-10T13:16:00.0" # http://odahub.io/ontology#EndTime | 33 T2 = "2004-12-20T13:16:00.0" # http://odahub.io/ontology#EndTime |
30 Radius = 2.5 # http://odahub.io/ontology#AngleDegrees | 34 Radius = 2.5 # http://odahub.io/ontology#AngleDegrees |
31 R_s = 0.2 # http://odahub.io/ontology#AngleDegrees | 35 R_s = 0.2 # http://odahub.io/ontology#AngleDegrees |
32 Emin = 100.0 # http://odahub.io/ontology#Energy_GeV | 36 Emin = 1 # http://odahub.io/ontology#Energy_TeV |
33 Emax = 10000.0 # http://odahub.io/ontology#Energy_GeV | 37 Emax = 100.0 # http://odahub.io/ontology#Energy_TeV |
34 NTbins = 10 # http://odahub.io/ontology#Integer | 38 NTbins = 30 # http://odahub.io/ontology#Integer |
35 | 39 |
36 _galaxy_wd = os.getcwd() | 40 _galaxy_wd = os.getcwd() |
37 | 41 |
38 with open("inputs.json", "r") as fd: | 42 with open("inputs.json", "r") as fd: |
39 inp_dic = json.load(fd) | 43 inp_dic = json.load(fd) |
103 flux = np.zeros(NTbins) | 107 flux = np.zeros(NTbins) |
104 flux_err = np.zeros(NTbins) | 108 flux_err = np.zeros(NTbins) |
105 flux_b = np.zeros(NTbins) | 109 flux_b = np.zeros(NTbins) |
106 flux_b_err = np.zeros(NTbins) | 110 flux_b_err = np.zeros(NTbins) |
107 Expos = np.zeros(NTbins) | 111 Expos = np.zeros(NTbins) |
112 src_cts = np.zeros(NTbins) | |
113 bkg_cts = np.zeros(NTbins) | |
108 for count, f in enumerate(OBSlist): | 114 for count, f in enumerate(OBSlist): |
115 pr.report_progress( | |
116 stage="Progress", progress=5.0 + 95 * count / len(OBSlist) | |
117 ) | |
109 hdul = fits.open("data/" + f) | 118 hdul = fits.open("data/" + f) |
110 RA_pnt = hdul[1].header["RA_PNT"] | 119 RA_pnt = hdul[1].header["RA_PNT"] |
111 DEC_pnt = hdul[1].header["DEC_PNT"] | 120 DEC_pnt = hdul[1].header["DEC_PNT"] |
112 Texp = hdul[1].header["LIVETIME"] | 121 Texp = hdul[1].header["LIVETIME"] |
113 Trun_start = hdul[1].header["TSTART"] | 122 Trun_start = hdul[1].header["TSTART"] |
122 ev = hdul["EVENTS"].data | 131 ev = hdul["EVENTS"].data |
123 ev_ra = ev["RA"] | 132 ev_ra = ev["RA"] |
124 ev_dec = ev["DEC"] | 133 ev_dec = ev["DEC"] |
125 ev_en = ev["ENERGY"] | 134 ev_en = ev["ENERGY"] |
126 ev_time = (ev["TIME"] - Trun_start) / 86400.0 + Tbegs[count] | 135 ev_time = (ev["TIME"] - Trun_start) / 86400.0 + Tbegs[count] |
136 Nevents = len(ev) | |
127 print(ev_time[0]) | 137 print(ev_time[0]) |
128 ev_coords = SkyCoord(ev_ra, ev_dec, unit="degree") | 138 ev_coords = SkyCoord(ev_ra, ev_dec, unit="degree") |
129 sep_s = ev_coords.separation(Coords_s).deg | 139 sep_s = ev_coords.separation(Coords_s).deg |
130 sep_b = ev_coords.separation(Coords_b).deg | 140 sep_b = ev_coords.separation(Coords_b).deg |
131 | 141 |
136 EEbins = np.concatenate((EEmin, [EEmax[-1]])) | 146 EEbins = np.concatenate((EEmin, [EEmax[-1]])) |
137 AA = hdu["EFFAREA"][0] + 1e-10 | 147 AA = hdu["EFFAREA"][0] + 1e-10 |
138 Thmin = hdu["THETA_LO"][0] | 148 Thmin = hdu["THETA_LO"][0] |
139 Thmax = hdu["THETA_HI"][0] | 149 Thmax = hdu["THETA_HI"][0] |
140 ind = np.argmin((Thmin - dist) ** 2) | 150 ind = np.argmin((Thmin - dist) ** 2) |
141 AA = AA[ind] * Texp * 1e4 | 151 mask = EE < Emin |
142 mask = np.where((sep_s < R_s)) | 152 ind_en = len(EE[mask]) |
143 cts1 = np.histogram2d(ev_time[mask], ev_en[mask], bins=[Tbins, EEbins])[0] | 153 Expos += np.histogram( |
144 mask = sep_b < R_s | 154 ev_time, |
145 cts2 = np.histogram2d(ev_time[mask], ev_en[mask], bins=[Tbins, EEbins])[0] | 155 weights=AA[ind, ind_en] * Texp * 1e4 * np.ones(Nevents) / Nevents, |
146 src = cts1 - cts2 | 156 bins=Tbins, |
147 src_err = sqrt(cts1 + cts2) | 157 )[0] |
148 flux += np.sum(src / AA, axis=1) | 158 mask = np.where((sep_s < R_s) & (ev_en > Emin) & (ev_en < Emax)) |
149 flux_err += np.sum(src_err / AA, axis=1) | 159 src_cts += np.histogram(ev_time[mask], bins=Tbins)[0] |
150 flux_b += np.sum(cts2 / AA, axis=1) | 160 mask = np.where((sep_b < R_s) & (ev_en > Emin) & (ev_en < Emax)) |
151 flux_b_err += np.sum(sqrt(cts2) / AA, axis=1) | 161 bkg_cts += np.histogram(ev_time[mask], bins=Tbins)[0] |
152 hdul.close() | 162 hdul.close() |
163 print(src_cts) | |
164 print(bkg_cts) | |
165 print(Expos) | |
166 src = src_cts - bkg_cts | |
167 src_err = sqrt(src_cts + bkg_cts) | |
168 flux = src / (Expos + 1) | |
169 flux_err = src_err / (Expos + 1) | |
170 flux_b = bkg_cts / (Expos + 1) | |
171 flux_b_err = sqrt(bkg_cts) / (Expos + 1) | |
172 # flux_err+=np.sum(src_err/AA,axis=1) | |
173 # flux_b+=np.sum(cts2/AA,axis=1) | |
174 # flux_b_err+=np.sum(sqrt(cts2)/AA,axis=1) | |
175 print(flux) | |
153 | 176 |
154 if message == "": | 177 if message == "": |
155 plt.errorbar( | 178 plt.errorbar( |
156 Tmean, | 179 Tmean, |
157 flux, | 180 flux, |
191 "Background[counts/cm2s]", | 214 "Background[counts/cm2s]", |
192 "Background_error[counts/cm2s]", | 215 "Background_error[counts/cm2s]", |
193 ) | 216 ) |
194 lc = ODAAstropyTable(Table(data, names=names)) | 217 lc = ODAAstropyTable(Table(data, names=names)) |
195 | 218 |
196 picture = bin_image # http://odahub.io/ontology#ODAPictureProduct | 219 png = bin_image # http://odahub.io/ontology#ODAPictureProduct |
197 lightcurve_astropy_table = lc # http://odahub.io/ontology#ODAAstropyTable | 220 table = lc # http://odahub.io/ontology#ODAAstropyTable |
198 | 221 |
199 # output gathering | 222 # output gathering |
200 _galaxy_meta_data = {} | 223 _galaxy_meta_data = {} |
201 _oda_outs = [] | 224 _oda_outs = [] |
202 _oda_outs.append(("out_Lightcurve_picture", "picture_galaxy.output", picture)) | 225 _oda_outs.append(("out_Lightcurve_png", "png_galaxy.output", png)) |
203 _oda_outs.append( | 226 _oda_outs.append(("out_Lightcurve_table", "table_galaxy.output", table)) |
204 ( | |
205 "out_Lightcurve_lightcurve_astropy_table", | |
206 "lightcurve_astropy_table_galaxy.output", | |
207 lightcurve_astropy_table, | |
208 ) | |
209 ) | |
210 | 227 |
211 for _outn, _outfn, _outv in _oda_outs: | 228 for _outn, _outfn, _outv in _oda_outs: |
212 _galaxy_outfile_name = os.path.join(_galaxy_wd, _outfn) | 229 _galaxy_outfile_name = os.path.join(_galaxy_wd, _outfn) |
213 if isinstance(_outv, str) and os.path.isfile(_outv): | 230 if isinstance(_outv, str) and os.path.isfile(_outv): |
214 shutil.move(_outv, _galaxy_outfile_name) | 231 shutil.move(_outv, _galaxy_outfile_name) |