Mercurial > repos > astroteam > hess_astro_tool

comparison Lightcurve.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=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)