Update tool AGN cacades to 0.0.1+galaxy0

tools/agn-cacades/.shed.yml

@@ -0,0 +1,9 @@
+categories:
+- Astronomy
+description: AGN cacades
+homepage_url: null
+long_description: AGN cacades
+name: agn_cacades_astro_tool
+owner: astroteam
+remote_repository_url: https://github.com/esg-epfl-apc/tools-astro/tree/main/tools
+type: unrestricted

tools/agn-cacades/Simulate.py

@@ -0,0 +1,314 @@
+#!/usr/bin/env python
+# coding: utf-8
+
+# flake8: noqa
+
+import json
+import os
+import shutil
+
+import matplotlib.pyplot as plt
+import numpy as np
+from numpy import log, pi
+from oda_api.data_products import ODAAstropyTable, PictureProduct
+from oda_api.json import CustomJSONEncoder
+from scipy.interpolate import interp1d
+
+get_ipython().system("pwd")   # noqa: F821
+
+get_ipython().system("git lfs install ; git lfs pull")   # noqa: F821
+
+get_ipython().run_cell_magic("bash", "", "./install.sh\n")   # noqa: F821
+
+Gamma = 2.6  # http://odahub.io/ontology#Float
+Emax = 1e15  # http://odahub.io/ontology#Float
+Ecut = 1e15  # http://odahub.io/ontology#Float
+B = 1e4  # http://odahub.io/ontology#Float
+source_size_cm = 3.0856e13  # http://odahub.io/ontology#Float
+background_norm_mode = "density_cm3"  # http://odahub.io/ontology#String ; oda:allowed_value "absolute","density_cm3","energy_density_eV_cm3"
+corona_backround_norm_cm3 = 4.18876e15  # http://odahub.io/ontology#Float
+disk_background_norm_cm3 = 1.39654e17  # http://odahub.io/ontology#Float
+min_steps = 100  # http://odahub.io/ontology#Integer ; oda:lower_limit 10 ; oda:upper_limit 10000
+
+_galaxy_wd = os.getcwd()
+
+with open("inputs.json", "r") as fd:
+    inp_dic = json.load(fd)
+if "_data_product" in inp_dic.keys():
+    inp_pdic = inp_dic["_data_product"]
+else:
+    inp_pdic = inp_dic
+
+for vn, vv in inp_pdic.items():
+    if vn != "_selector":
+        globals()[vn] = type(globals()[vn])(vv)
+
+background_norm_mode_index = [
+    "absolute",
+    "density_cm3",
+    "energy_density_eV_cm3",
+].index(background_norm_mode)
+distance_Mpc = source_size_cm / 3.0856775807e24
+electron_distance_Mpc = 100 * distance_Mpc  # electrons propagate longer
+min_step_Mpc = distance_Mpc / min_steps
+electron_min_step_Mpc = electron_distance_Mpc / min_steps
+
+# fname='pp_a'+str(Gamma)+'_'+str(f"{Ecut:.2e}")+'_'+str(B)
+fname = "run"
+get_ipython().system(   # noqa: F821
+    "cp propagation/pp_a2.4_E1e14.xsw propagation/{fname}.xsw"
+)
+
+get_ipython().system(   # noqa: F821
+    "python propagation/replacePar CommonAlpha {str(Gamma)} propagation/{fname}.xsw"
+)
+get_ipython().system(   # noqa: F821
+    "python propagation/getPar CommonAlpha propagation/{fname}.xsw"
+)
+
+get_ipython().system(   # noqa: F821
+    "python propagation/replacePar Emax {str(Emax)} propagation/{fname}.xsw"
+)
+get_ipython().system("python propagation/getPar Emax propagation/{fname}.xsw")   # noqa: F821
+
+Ecut_rel = Ecut / Emax
+get_ipython().system(   # noqa: F821
+    "python propagation/replacePar injSpectraHigherCutoff {str(Ecut_rel)} propagation/{fname}.xsw"
+)
+get_ipython().system(   # noqa: F821
+    "python propagation/getPar injSpectraHigherCutoff propagation/{fname}.xsw"
+)
+
+get_ipython().system(   # noqa: F821
+    "python propagation/replacePar B0 {str(B)} propagation/step??_template.xsw"
+)
+get_ipython().system(   # noqa: F821
+    "python propagation/getPar B0 propagation/step2g_template.xsw"
+)
+
+get_ipython().system(   # noqa: F821
+    "python propagation/replacePar iroMult {str(disk_background_norm_cm3)} propagation/step??_template.xsw"
+)
+get_ipython().system(   # noqa: F821
+    "python propagation/getPar iroMult propagation/step2g_template.xsw"
+)
+
+get_ipython().system(   # noqa: F821
+    "python propagation/replacePar IROBackgroundNormMode {str(background_norm_mode_index)} propagation/step??_template.xsw"
+)
+get_ipython().system(   # noqa: F821
+    "python propagation/getPar IROBackgroundNormMode propagation/step2g_template.xsw"
+)
+
+get_ipython().system(   # noqa: F821
+    "python propagation/replacePar CustomBackgroundNormMode {str(background_norm_mode_index)} propagation/step??_template.xsw"
+)
+get_ipython().system(   # noqa: F821
+    "python propagation/getPar CustomBackgroundNormMode propagation/step2g_template.xsw"
+)
+
+get_ipython().system(   # noqa: F821
+    "python propagation/replacePar L_short_distance_test {str(distance_Mpc)} propagation/step?g_template.xsw"
+)
+get_ipython().system(   # noqa: F821
+    "python propagation/getPar L_short_distance_test propagation/step2g_template.xsw"
+)
+
+get_ipython().system(   # noqa: F821
+    "python propagation/replacePar L_short_distance_test {str(electron_distance_Mpc)} propagation/step?e_template.xsw"
+)
+get_ipython().system(   # noqa: F821
+    "python propagation/getPar L_short_distance_test propagation/step3e_template.xsw"
+)
+
+get_ipython().system(   # noqa: F821
+    "python propagation/replacePar microStep {str(min_step_Mpc)} propagation/step?g_template.xsw"
+)
+get_ipython().system(   # noqa: F821
+    "python propagation/getPar microStep propagation/step2g_template.xsw"
+)
+
+get_ipython().system(   # noqa: F821
+    "python propagation/replacePar microStep {str(electron_min_step_Mpc)} propagation/step?e_template.xsw"
+)
+get_ipython().system(   # noqa: F821
+    "python propagation/getPar microStep propagation/step3e_template.xsw"
+)
+
+get_ipython().run_cell_magic(   # noqa: F821
+    "bash", "", "cd propagation\n./run_all.sh run.xsw\n"
+)
+
+get_ipython().system("cat propagation/results/run_final/tot")   # noqa: F821
+
+d = np.genfromtxt("propagation/results/run_final/tot")
+E = d[:, 0]
+gam = d[:, 1]
+nu = d[:, 2]
+elec = d[:, 3]
+prot = d[:, 4]
+
+factor = 4
+plt.plot(E / 1e9, gam * factor)
+plt.plot(E / 1e9, elec * factor)
+plt.plot(E / 1e9, nu * factor)
+
+d = np.genfromtxt("NGC_1068_contour.csv")
+gam = d[:, 0]
+f = d[:, 1]
+x_nu = np.logspace(3, 5, 10)
+ymax = np.zeros(len(x_nu))
+ymin = np.ones(len(x_nu))
+for i in range(len(gam)):
+    y = 3 * f[i] * 1e-11 * (x_nu / 1e3) ** (2 - gam[i])
+    ymax = np.maximum(y, ymax)
+    ymin = np.minimum(y, ymin)
+# plt.fill_between(x_nu,ymin,ymax,color='blue',alpha=0.5)
+
+d = np.genfromtxt("SED_1068.csv")
+nu = d[:, 0]
+nufnu = d[:, 1]
+en = 2 * pi * 6.6e-16 * nu / 1e9
+efe = nufnu * 1e-23 / 1.6
+# plt.scatter(en,efe)
+
+d = np.genfromtxt("Fermi_1068.csv")
+ee = d[:, 0]
+ff = d[:, 1]
+ff_max = d[:, 2]
+ff_min = d[:, 3]
+# plt.errorbar(ee,ff,yerr=[ff-ff_min,ff_max-ff])
+
+A = 7e43 / 1.6e-12 / 1e6 / log(10 / 2.0)  # 1/eV/s
+FX0 = A * 1e6 / 1e12 / 4 / pi / (16.3 * 3e24) ** 2
+print(FX0)
+x = [1e-6, 1e-4]
+y = [FX0, FX0]
+# plt.plot(x,y,color='black',linestyle='dashed')
+A = 4e43 / 1.6e-12 / 1e6 / log(10 / 2.0)  # 1/eV/s
+FX0 = A * 1e6 / 1e12 / 4 / pi / (16.3 * 3e24) ** 2
+A = 14e43 / 1.6e-12 / 1e6 / log(10 / 2.0)  # 1/eV/s
+FX1 = A * 1e6 / 1e12 / 4 / pi / (16.3 * 3e24) ** 2
+x = [1e-6, 1e-4, 1e-4, 1e-6]
+y = [FX1, FX1, FX0, FX0]
+# plt.fill(x,y,alpha=0.3,color='black',linestyle='none')
+
+plt.xscale("log")
+# plt.yscale('log')
+plt.ylim(1e-14, 1e-8)
+plt.xlim(1e-12, 1e6)
+
+d = np.genfromtxt("NGC_1068_contour.csv")
+gam = d[:, 0]
+f = d[:, 1]
+x_nu = np.logspace(3, 5, 10)
+ymax = np.zeros(len(x_nu))
+ymin = np.ones(len(x_nu))
+for i in range(len(gam)):
+    y = 3 * f[i] * 1e-11 * (x_nu / 1e3) ** (2 - gam[i])
+    ymax = np.maximum(y, ymax)
+    ymin = np.minimum(y, ymin)
+plt.fill_between(x_nu, ymin, ymax, color="blue", alpha=0.5)
+
+d = np.genfromtxt("SED_1068.csv")
+nu = d[:, 0]
+nufnu = d[:, 1]
+en = 2 * pi * 6.6e-16 * nu / 1e9
+efe = nufnu * 1e-23 / 1.6
+plt.scatter(en, efe)
+
+d = np.genfromtxt("Fermi_1068.csv")
+ee = d[:, 0]
+ff = d[:, 1]
+ff_max = d[:, 2]
+ff_min = d[:, 3]
+plt.errorbar(ee, ff, yerr=[ff - ff_min, ff_max - ff])
+
+A = 7e43 / 1.6e-12 / 1e6 / log(10 / 2.0)  # 1/eV/s
+FX0 = A * 1e6 / 1e12 / 4 / pi / (16.3 * 3e24) ** 2
+print(FX0)
+x = [1e-6, 1e-4]
+y = [FX0, FX0]
+plt.plot(x, y, color="black", linestyle="dashed")
+A = 4e43 / 1.6e-12 / 1e6 / log(10 / 2.0)  # 1/eV/s
+FX0 = A * 1e6 / 1e12 / 4 / pi / (16.3 * 3e24) ** 2
+A = 14e43 / 1.6e-12 / 1e6 / log(10 / 2.0)  # 1/eV/s
+FX1 = A * 1e6 / 1e12 / 4 / pi / (16.3 * 3e24) ** 2
+x = [1e-6, 1e-4, 1e-4, 1e-6]
+y = [FX1, FX1, FX0, FX0]
+plt.fill(x, y, alpha=0.3, color="black", linestyle="none")
+
+spec_file = "propagation/results/run_final/tot"
+d = np.genfromtxt(spec_file)
+E = d[:, 0] * 1e-9
+gam = d[:, 1]  # for original format d[:,3]
+nu = d[:, 2]  # for original format (d[:,4]+d[:,5]+d[:,7]+d[:,8])
+f_nu = interp1d(np.log(E), nu, kind="linear", fill_value=0.0)
+interp_nu = f_nu(np.log(x_nu))
+N = 0.8 * np.min(ymax / interp_nu)
+plt.plot(E, N * gam * (E > 5.0e-9), label="$\\gamma$", color="blue", alpha=0.3)
+plt.plot(E, N * gam * (E > 5.0e-7), color="black")
+plt.plot(E, N * nu, label="$\\nu$")
+backgr = np.genfromtxt("propagation/results/run_step4g/fin/backgr")
+Eb = backgr[:, 0] * 1e-9
+FE2 = (
+    backgr[:, 0]
+    * backgr[:, 1]
+    * 1e-12
+    * 3e10
+    / 4
+    / np.pi
+    * (1e13 / 4.6285163711e25) ** 2
+)
+plt.plot(Eb, 20 * FE2, label="backgr3")
+
+plt.xscale("log")
+plt.yscale("log")
+plt.ylim(6e-14, 3e-8)
+plt.xlim(3e-11, 1e6)
+plt.tick_params(labelsize=12)
+plt.xlabel(r"$E$, GeV", fontsize=12)
+plt.ylabel(r"$EF_E$, TeV/(cm$^2$s)", fontsize=12)
+plt.savefig("Spectrum.png", format="png", bbox_inches="tight")
+
+bin_image = PictureProduct.from_file("Spectrum.png")
+from astropy.table import Table
+
+data = [E, gam, nu]
+names = ("E[eV]", "E2 dN/dE gamma [TeV/cm2s]", "E2 dN/dE nu[ TeV/cm2 s]")
+spectrum = ODAAstropyTable(Table(data, names=names))
+
+picture = bin_image  # http://odahub.io/ontology#ODAPictureProduct
+spectrum_astropy_table = spectrum  # http://odahub.io/ontology#ODAAstropyTable
+
+# output gathering
+_galaxy_meta_data = {}
+_oda_outs = []
+_oda_outs.append(("out_Simulate_picture", "picture_galaxy.output", picture))
+_oda_outs.append(
+    (
+        "out_Simulate_spectrum_astropy_table",
+        "spectrum_astropy_table_galaxy.output",
+        spectrum_astropy_table,
+    )
+)
+
+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 ***")

tools/agn-cacades/agn_cacades_astro_tool.xml

@@ -0,0 +1,63 @@
+<tool id="agn_cacades_astro_tool" name="AGN cacades" version="0.0.1+galaxy0" profile="23.0">
+  <requirements>
+    <requirement type="package" version="3.29.6">cmake</requirement>
+    <requirement type="package" version="12.3.0">gfortran</requirement>
+    <requirement type="package" version="2.7">gsl</requirement>
+    <requirement type="package" version="1.7.0">compilers</requirement>
+    <requirement type="package" version="3.2.5">xerces-c</requirement>
+    <requirement type="package" version="2.2.2">pandas</requirement>
+    <requirement type="package" version="1.2.20">oda-api</requirement>
+    <requirement type="package" version="3.9.0">matplotlib</requirement>
+    <requirement type="package" version="2.32.3">requests</requirement>
+    <requirement type="package" version="0.4.7">astroquery</requirement>
+    <requirement type="package" version="7.0.0">rdflib</requirement>
+    <requirement type="package" version="4.12.3">beautifulsoup4</requirement>
+    <requirement type="package" version="2.0.0">sparqlwrapper</requirement>
+    <requirement type="package" version="3.3">networkx</requirement>
+    <requirement type="package" version="1.14.0">scipy</requirement>
+    <requirement type="package" version="8.26.0">ipython</requirement>
+    <requirement type="package" version="6.1.1">astropy</requirement>
+    <!--Requirements string 'nb2workflow[cwl,service,rdf,mmoda]>=1.3.30 
+' can't be converted automatically. Please add the galaxy/conda requirement manually or modify the requirements file!-->
+    <requirement type="package" version="7.16.4">nbconvert</requirement>
+  </requirements>
+  <command detect_errors="exit_code">ipython '$__tool_directory__/Simulate.py'</command>
+  <configfiles>
+    <inputs name="inputs" filename="inputs.json" data_style="paths" />
+  </configfiles>
+  <inputs>
+    <param name="Gamma" type="float" value="2.6" label="Gamma" />
+    <param name="Emax" type="float" value="1000000000000000.0" label="Emax" />
+    <param name="Ecut" type="float" value="1000000000000000.0" label="Ecut" />
+    <param name="B" type="float" value="10000.0" label="B" />
+    <param name="source_size_cm" type="float" value="30856000000000.0" label="source_size_cm" />
+    <param name="background_norm_mode" type="select" label="background_norm_mode">
+      <option value="absolute">absolute</option>
+      <option value="density_cm3" selected="true">density_cm3</option>
+      <option value="energy_density_eV_cm3">energy_density_eV_cm3</option>
+    </param>
+    <param name="corona_backround_norm_cm3" type="float" value="4188760000000000.0" label="corona_backround_norm_cm3" />
+    <param name="disk_background_norm_cm3" type="float" value="1.39654e+17" label="disk_background_norm_cm3" />
+    <param name="min_steps" type="integer" value="100" label="min_steps" min="10" max="10000" />
+  </inputs>
+  <outputs>
+    <data label="${tool.name} -&gt; Simulate picture" name="out_Simulate_picture" format="auto" from_work_dir="picture_galaxy.output" />
+    <data label="${tool.name} -&gt; Simulate spectrum_astropy_table" name="out_Simulate_spectrum_astropy_table" format="auto" from_work_dir="spectrum_astropy_table_galaxy.output" />
+  </outputs>
+  <tests>
+    <test expect_num_outputs="2">
+      <param name="Gamma" value="2.6" />
+      <param name="Emax" value="1000000000000000.0" />
+      <param name="Ecut" value="1000000000000000.0" />
+      <param name="B" value="10000.0" />
+      <param name="source_size_cm" value="30856000000000.0" />
+      <param name="background_norm_mode" value="density_cm3" />
+      <param name="corona_backround_norm_cm3" value="4188760000000000.0" />
+      <param name="disk_background_norm_cm3" value="1.39654e+17" />
+      <param name="min_steps" value="100" />
+      <assert_stdout>
+        <has_text text="*** Job finished successfully ***" />
+      </assert_stdout>
+    </test>
+  </tests>
+</tool>