Mercurial > repos > astroteam > crbeam_astro_tool
diff Generate_events.py @ 0:f40d05521dca draft default tip
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 |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Generate_events.py Fri Apr 25 19:33:20 2025 +0000 @@ -0,0 +1,249 @@ +#!/usr/bin/env python +# coding: utf-8 + +#!/usr/bin/env python + +# This script is generated with nb2galaxy + +# flake8: noqa + +import json +import os +import shutil + +from oda_api.json import CustomJSONEncoder + +# oda:usesResource oda:CRBeamS3 . +# oda:CRBeamS3 a oda:S3 . +# oda:CRBeamS3 oda:resourceBindingEnvVarName "S3_CREDENTIALS" . + +src_name = "1ES 1215+303" # http://odahub.io/ontology#AstrophysicalObject + +z_start = 0.13 # http://odahub.io/ontology#Float +Npart = 10000 # http://odahub.io/ontology#Integer ; oda:lower_limit 1 ; oda:upper_limit 100000 +particle_type = "gamma" # http://odahub.io/ontology#String ; oda:allowed_value "gamma","electron","proton" +Emax = 50 # http://odahub.io/ontology#Energy_TeV +Emin = 0.01 # http://odahub.io/ontology#Energy_TeV +EminSource = 0.01 # http://odahub.io/ontology#Energy_TeV +Gamma = 2.0 # http://odahub.io/ontology#Float +EGMF_fG = 100 # http://odahub.io/ontology#Float +lmaxEGMF_Mpc = 5 # http://odahub.io/ontology#Float +jet_half_size = 180.0 # http://odahub.io/ontology#AngleDegrees +jet_direction = 0.0 # http://odahub.io/ontology#AngleDegrees +psf = 180.0 # http://odahub.io/ontology#AngleDegrees +window_size_RA = 4.0 # http://odahub.io/ontology#AngleDegrees +window_size_DEC = 4.0 # http://odahub.io/ontology#AngleDegrees +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" + +_galaxy_wd = os.getcwd() + +with open("inputs.json", "r") as fd: + inp_dic = json.load(fd) +if "C_data_product_" in inp_dic.keys(): + inp_pdic = inp_dic["C_data_product_"] +else: + inp_pdic = inp_dic +src_name = str(inp_pdic["src_name"]) +z_start = float(inp_pdic["z_start"]) +Npart = int(inp_pdic["Npart"]) +particle_type = str(inp_pdic["particle_type"]) +Emax = float(inp_pdic["Emax"]) +Emin = float(inp_pdic["Emin"]) +EminSource = float(inp_pdic["EminSource"]) +Gamma = float(inp_pdic["Gamma"]) +EGMF_fG = float(inp_pdic["EGMF_fG"]) +lmaxEGMF_Mpc = float(inp_pdic["lmaxEGMF_Mpc"]) +jet_half_size = float(inp_pdic["jet_half_size"]) +jet_direction = float(inp_pdic["jet_direction"]) +psf = float(inp_pdic["psf"]) +window_size_RA = float(inp_pdic["window_size_RA"]) +window_size_DEC = float(inp_pdic["window_size_DEC"]) +EBL = str(inp_pdic["EBL"]) + +import os + +from astropy.coordinates import SkyCoord +from astropy.table import Table +from oda_api.api import ProgressReporter +from oda_api.data_products import ODAAstropyTable +from utils import find_redshift + +if z_start <= 0: + z_start = find_redshift(src_name) + +source = SkyCoord.from_name(src_name, frame="icrs", parse=False, cache=True) + +pr = ProgressReporter() +pr.report_progress( + stage="Initialization", + progress=0, + substage="spectra", + subprogress=30, + message="some message", +) + +import light_curve as lc +import numpy as np +from crbeam import CRbeam + +EGMF = EGMF_fG * 1e-15 + +# internal units are eV +Emax *= 1e12 +Emin *= 1e12 +EminSource *= 1e12 + +background = { + "Franceschini 2017": 12, + "Stecker 2016 lower limit": 10, + "Stecker 2016 upper limit": 11, + "Inoue 2012 Baseline": 3, + "Inoue 2012 lower limit": 4, + "Inoue 2012 upper limit": 5, +}[EBL] + +prog = CRbeam( + z=z_start, + nparticles=Npart, + primary=particle_type, + emax=Emax, + emin=Emin, + emin_source=EminSource, + EGMF=EGMF, + lmaxEGMF=lmaxEGMF_Mpc, + background=background, +) +cmd = prog.command +cmd + +# Here is one way to call CRbeam without global cache +# prog.run(force_overwrite=False) +# Here we call CRbeam with global cache +# prog.run_cached(overwrite_local_cache=True) +# to clear s3 cache run the following command +# prog.remove_cache() + +# To report the progress we will split running CRbeam into 10 parts +n_steps = 10 +# Initialize multistep simulation +data_exists = not prog.start_multistage_run( + overwrite_local_cache=True, n_steps=n_steps +) +proceed = not data_exists + +if proceed: + for step in range(n_steps): + pr.report_progress( + stage="Running Simulation", progress=100 * step // n_steps + ) + proceed = prog.simulation_step() + # todo: report progress using rest API + pr.report_progress(stage="Running Simulation", progress=100) + +assert not proceed, "must be completed before this cell" +if not data_exists: + prog.end_multistep_run() + +def adjust_weights(mc_file, power): + # converting weights to mimic required injection spectrum power + header = "" + with open(mc_file, "rt") as lines: + for line in lines: + if len(line) > 0 and line[0] == "#": + header += line[1:].strip() + "\n" + else: + break + weight_col = 2 + E_src_col = 12 + data = np.loadtxt(mc_file) + weight = data[:, weight_col - 1] + E_src = data[:, E_src_col - 1] + orig_power = ( + prog.power_law + ) # CRBeam is always called with fixed power=1 to optimize cache usage + weight *= np.power(E_src / Emax, -(power - orig_power)) + output_file = f"{mc_file}_p{power}" + np.savetxt(output_file, data, header=header.strip(), fmt="%.6g") + return output_file + +# this code will not execute program since files already exist on server +# prog.run(force_overwrite=False) + +# one can upload cache explicitely by call +# prog.upload_cache() + +pr.report_progress(stage="Building Plots and Tables", progress=0) + +print(prog.output_path) +get_ipython().system("ls {prog.output_path}") # noqa: F821 + +mc_file = prog.output_path + "/photon" + +if Emax != EminSource: + mc_file = adjust_weights(mc_file, Gamma) + +# rotating the beam + +if EGMF > 0: + from mc_rotate import mc_rotate + + mc_rotated_file = mc_rotate( + mc_file, jet_half_size, jet_direction, psf_deg=psf + ) +else: + mc_rotated_file = mc_file +mc_rotated_file + +# reading the source distance in Mpc from the data file +T_Mpc = lc.get_distance_Mpc(mc_file) +T_Mpc + +from oda_api.data_products import ODAAstropyTable + +d = np.genfromtxt(mc_rotated_file, skip_header=3) +# E/eV,weight,Theta,Phi,dT_raw/T dT_calc/T,z_cascade_production,N_interactions,z_src,E_src/eV,z_prod +names = ( + "E/eV", + "weight", + "Theta", + "Phi", + "dT_raw/T", + "dT_calc/T", + "z_cascade_production", + "N_interactions", + "z_src", + "E_src/eV", + "z_prod", +) +res = ODAAstropyTable(Table(d, names=names)) + +pr.report_progress(stage="Building Plots and Tables", progress=100) + +Event_file = res # http://odahub.io/ontology#ODAAstropyTable + +# output gathering +_galaxy_meta_data = {} +_oda_outs = [] +_oda_outs.append( + ("out_Generate_events_Event_file", "Event_file_galaxy.output", Event_file) +) + +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 ***")