Mercurial > repos > astroteam > crbeam_astro_tool

comparison Generate_events.py @ 0:f40d05521dca 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 de01e3c02a26cd6353a6b9b6f8d1be44de8ccd54
author astroteam
date Fri, 25 Apr 2025 19:33:20 +0000
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:f40d05521dca
1 #!/usr/bin/env python
2 # coding: utf-8
3
4 #!/usr/bin/env python
5
6 # This script is generated with nb2galaxy
7
8 # flake8: noqa
9
10 import json
11 import os
12 import shutil
13
14 from oda_api.json import CustomJSONEncoder
15
16 # oda:usesResource oda:CRBeamS3 .
17 # oda:CRBeamS3 a oda:S3 .
18 # oda:CRBeamS3 oda:resourceBindingEnvVarName "S3_CREDENTIALS" .
19
20 src_name = "1ES 1215+303" # http://odahub.io/ontology#AstrophysicalObject
21
22 z_start = 0.13 # http://odahub.io/ontology#Float
23 Npart = 10000 # http://odahub.io/ontology#Integer ; oda:lower_limit 1 ; oda:upper_limit 100000
24 particle_type = "gamma" # http://odahub.io/ontology#String ; oda:allowed_value "gamma","electron","proton"
25 Emax = 50 # http://odahub.io/ontology#Energy_TeV
26 Emin = 0.01 # http://odahub.io/ontology#Energy_TeV
27 EminSource = 0.01 # http://odahub.io/ontology#Energy_TeV
28 Gamma = 2.0 # http://odahub.io/ontology#Float
29 EGMF_fG = 100 # http://odahub.io/ontology#Float
30 lmaxEGMF_Mpc = 5 # http://odahub.io/ontology#Float
31 jet_half_size = 180.0 # http://odahub.io/ontology#AngleDegrees
32 jet_direction = 0.0 # http://odahub.io/ontology#AngleDegrees
33 psf = 180.0 # http://odahub.io/ontology#AngleDegrees
34 window_size_RA = 4.0 # http://odahub.io/ontology#AngleDegrees
35 window_size_DEC = 4.0 # http://odahub.io/ontology#AngleDegrees
36 EBL = "Franceschini 2017" # http://odahub.io/ontology#String ; oda:allowed_value "Franceschini 2017","Stecker 2016 lower limit","Stecker 2016 upper limit","Inoue 2012 Baseline","Inoue 2012 lower limit","Inoue 2012 upper limit"
37
38 _galaxy_wd = os.getcwd()
39
40 with open("inputs.json", "r") as fd:
41 inp_dic = json.load(fd)
42 if "C_data_product_" in inp_dic.keys():
43 inp_pdic = inp_dic["C_data_product_"]
44 else:
45 inp_pdic = inp_dic
46 src_name = str(inp_pdic["src_name"])
47 z_start = float(inp_pdic["z_start"])
48 Npart = int(inp_pdic["Npart"])
49 particle_type = str(inp_pdic["particle_type"])
50 Emax = float(inp_pdic["Emax"])
51 Emin = float(inp_pdic["Emin"])
52 EminSource = float(inp_pdic["EminSource"])
53 Gamma = float(inp_pdic["Gamma"])
54 EGMF_fG = float(inp_pdic["EGMF_fG"])
55 lmaxEGMF_Mpc = float(inp_pdic["lmaxEGMF_Mpc"])
56 jet_half_size = float(inp_pdic["jet_half_size"])
57 jet_direction = float(inp_pdic["jet_direction"])
58 psf = float(inp_pdic["psf"])
59 window_size_RA = float(inp_pdic["window_size_RA"])
60 window_size_DEC = float(inp_pdic["window_size_DEC"])
61 EBL = str(inp_pdic["EBL"])
62
63 import os
64
65 from astropy.coordinates import SkyCoord
66 from astropy.table import Table
67 from oda_api.api import ProgressReporter
68 from oda_api.data_products import ODAAstropyTable
69 from utils import find_redshift
70
71 if z_start <= 0:
72 z_start = find_redshift(src_name)
73
74 source = SkyCoord.from_name(src_name, frame="icrs", parse=False, cache=True)
75
76 pr = ProgressReporter()
77 pr.report_progress(
78 stage="Initialization",
79 progress=0,
80 substage="spectra",
81 subprogress=30,
82 message="some message",
83 )
84
85 import light_curve as lc
86 import numpy as np
87 from crbeam import CRbeam
88
89 EGMF = EGMF_fG * 1e-15
90
91 # internal units are eV
92 Emax *= 1e12
93 Emin *= 1e12
94 EminSource *= 1e12
95
96 background = {
97 "Franceschini 2017": 12,
98 "Stecker 2016 lower limit": 10,
99 "Stecker 2016 upper limit": 11,
100 "Inoue 2012 Baseline": 3,
101 "Inoue 2012 lower limit": 4,
102 "Inoue 2012 upper limit": 5,
103 }[EBL]
104
105 prog = CRbeam(
106 z=z_start,
107 nparticles=Npart,
108 primary=particle_type,
109 emax=Emax,
110 emin=Emin,
111 emin_source=EminSource,
112 EGMF=EGMF,
113 lmaxEGMF=lmaxEGMF_Mpc,
114 background=background,
115 )
116 cmd = prog.command
117 cmd
118
119 # Here is one way to call CRbeam without global cache
120 # prog.run(force_overwrite=False)
121 # Here we call CRbeam with global cache
122 # prog.run_cached(overwrite_local_cache=True)
123 # to clear s3 cache run the following command
124 # prog.remove_cache()
125
126 # To report the progress we will split running CRbeam into 10 parts
127 n_steps = 10
128 # Initialize multistep simulation
129 data_exists = not prog.start_multistage_run(
130 overwrite_local_cache=True, n_steps=n_steps
131 )
132 proceed = not data_exists
133
134 if proceed:
135 for step in range(n_steps):
136 pr.report_progress(
137 stage="Running Simulation", progress=100 * step // n_steps
138 )
139 proceed = prog.simulation_step()
140 # todo: report progress using rest API
141 pr.report_progress(stage="Running Simulation", progress=100)
142
143 assert not proceed, "must be completed before this cell"
144 if not data_exists:
145 prog.end_multistep_run()
146
147 def adjust_weights(mc_file, power):
148 # converting weights to mimic required injection spectrum power
149 header = ""
150 with open(mc_file, "rt") as lines:
151 for line in lines:
152 if len(line) > 0 and line[0] == "#":
153 header += line[1:].strip() + "\n"
154 else:
155 break
156 weight_col = 2
157 E_src_col = 12
158 data = np.loadtxt(mc_file)
159 weight = data[:, weight_col - 1]
160 E_src = data[:, E_src_col - 1]
161 orig_power = (
162 prog.power_law
163 ) # CRBeam is always called with fixed power=1 to optimize cache usage
164 weight *= np.power(E_src / Emax, -(power - orig_power))
165 output_file = f"{mc_file}_p{power}"
166 np.savetxt(output_file, data, header=header.strip(), fmt="%.6g")
167 return output_file
168
169 # this code will not execute program since files already exist on server
170 # prog.run(force_overwrite=False)
171
172 # one can upload cache explicitely by call
173 # prog.upload_cache()
174
175 pr.report_progress(stage="Building Plots and Tables", progress=0)
176
177 print(prog.output_path)
178 get_ipython().system("ls {prog.output_path}") # noqa: F821
179
180 mc_file = prog.output_path + "/photon"
181
182 if Emax != EminSource:
183 mc_file = adjust_weights(mc_file, Gamma)
184
185 # rotating the beam
186
187 if EGMF > 0:
188 from mc_rotate import mc_rotate
189
190 mc_rotated_file = mc_rotate(
191 mc_file, jet_half_size, jet_direction, psf_deg=psf
192 )
193 else:
194 mc_rotated_file = mc_file
195 mc_rotated_file
196
197 # reading the source distance in Mpc from the data file
198 T_Mpc = lc.get_distance_Mpc(mc_file)
199 T_Mpc
200
201 from oda_api.data_products import ODAAstropyTable
202
203 d = np.genfromtxt(mc_rotated_file, skip_header=3)
204 # E/eV,weight,Theta,Phi,dT_raw/T dT_calc/T,z_cascade_production,N_interactions,z_src,E_src/eV,z_prod
205 names = (
206 "E/eV",
207 "weight",
208 "Theta",
209 "Phi",
210 "dT_raw/T",
211 "dT_calc/T",
212 "z_cascade_production",
213 "N_interactions",
214 "z_src",
215 "E_src/eV",
216 "z_prod",
217 )
218 res = ODAAstropyTable(Table(d, names=names))
219
220 pr.report_progress(stage="Building Plots and Tables", progress=100)
221
222 Event_file = res # http://odahub.io/ontology#ODAAstropyTable
223
224 # output gathering
225 _galaxy_meta_data = {}
226 _oda_outs = []
227 _oda_outs.append(
228 ("out_Generate_events_Event_file", "Event_file_galaxy.output", Event_file)
229 )
230
231 for _outn, _outfn, _outv in _oda_outs:
232 _galaxy_outfile_name = os.path.join(_galaxy_wd, _outfn)
233 if isinstance(_outv, str) and os.path.isfile(_outv):
234 shutil.move(_outv, _galaxy_outfile_name)
235 _galaxy_meta_data[_outn] = {"ext": "_sniff_"}
236 elif getattr(_outv, "write_fits_file", None):
237 _outv.write_fits_file(_galaxy_outfile_name)
238 _galaxy_meta_data[_outn] = {"ext": "fits"}
239 elif getattr(_outv, "write_file", None):
240 _outv.write_file(_galaxy_outfile_name)
241 _galaxy_meta_data[_outn] = {"ext": "_sniff_"}
242 else:
243 with open(_galaxy_outfile_name, "w") as fd:
244 json.dump(_outv, fd, cls=CustomJSONEncoder)
245 _galaxy_meta_data[_outn] = {"ext": "json"}
246
247 with open(os.path.join(_galaxy_wd, "galaxy.json"), "w") as fd:
248 json.dump(_galaxy_meta_data, fd)
249 print("*** Job finished successfully ***")