From 632117266ccaac36b84fe699d8e40cc5bc69334d Mon Sep 17 00:00:00 2001
From: ODA bot <noreply@odahub.io>
Date: Wed, 17 Jan 2024 15:14:53 +0000
Subject: [PATCH] update tool HESS

---
 tools/hess/Image.py            | 252 +++++++++++++++++++++++++++++
 tools/hess/Lightcurve.py       | 280 +++++++++++++++++++++++++++++++++
 tools/hess/Spectrum.py         | 257 ++++++++++++++++++++++++++++++
 tools/hess/hess_astro_tool.xml | 138 ++++++++++++++++
 4 files changed, 927 insertions(+)
 create mode 100644 tools/hess/Image.py
 create mode 100644 tools/hess/Lightcurve.py
 create mode 100644 tools/hess/Spectrum.py
 create mode 100644 tools/hess/hess_astro_tool.xml

diff --git a/tools/hess/Image.py b/tools/hess/Image.py
new file mode 100644
index 00000000..ac108f0b
--- /dev/null
+++ b/tools/hess/Image.py
@@ -0,0 +1,252 @@
+#!/usr/bin/env python
+# coding: utf-8
+
+# In[ ]:
+
+
+import json
+import os
+import shutil
+import sys
+
+try:
+    import numpy as np
+    _numpy_available = True
+except ImportError:
+    _numpy_available = False
+
+try:
+    from oda_api.json import CustomJSONEncoder
+except ImportError:
+    from json import JSONEncoder as CustomJSONEncoder
+
+_galaxy_wd = os.getcwd()
+
+
+# In[1]:
+
+
+import astropy.units as u
+from astropy.coordinates import SkyCoord
+
+import matplotlib.pyplot as plt
+from gammapy.data import DataStore
+from gammapy.makers import MapDatasetMaker
+from gammapy.makers.utils import make_theta_squared_table
+from gammapy.maps import Map, MapAxis, WcsGeom
+from astropy.io import fits
+import numpy as np
+from numpy import pi,cos,sin,sqrt
+import os
+from astropy import wcs
+from astropy.io import fits
+from oda_api.data_products import PictureProduct
+from oda_api.data_products import ImageDataProduct
+from astropy.time import Time
+
+
+# In[2]:
+
+
+if(os.path.exists('hess_dl3_dr1.tar.gz')==False):
+    get_ipython().system('wget https://zenodo.org/record/1421099/files/hess_dl3_dr1.tar.gz')
+    get_ipython().system('tar -zxvf hess_dl3_dr1.tar.gz')
+
+
+# In[3]:
+
+
+src_name='Crab' #http://odahub.io/ontology#AstrophysicalObject
+RA=83.628700  # http://odahub.io/ontology#PointOfInterestRA
+DEC = 22.014700 # http://odahub.io/ontology#PointOfInterestDEC
+T1='2000-10-09T13:16:00.0'# http://odahub.io/ontology#StartTime
+T2='2022-10-10T13:16:00.0' # http://odahub.io/ontology#EndTime
+Radius=2.5  #http://odahub.io/ontology#AngleDegrees
+pixsize=0.1 #http://odahub.io/ontology#AngleDegrees
+Emin=100.    #http://odahub.io/ontology#Energy_GeV
+Emax=10000. #http://odahub.io/ontology#Energy_GeV
+
+
+# In[ ]:
+
+
+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)
+
+
+# In[4]:
+
+
+T1=Time(T1, format='isot', scale='utc').mjd
+T2=Time(T2, format='isot', scale='utc').mjd
+message=''
+RA_pnts=[]
+DEC_pnts=[]
+DL3_files=[]
+OBSIDs=[]
+Tstart=[]
+Tstop=[]
+flist=os.listdir('data')
+for f in flist:
+    if(f[-7:]=='fits.gz'):
+        DL3_files.append(f)
+        OBSIDs.append(int(f[20:26]))
+        hdul=fits.open('data/'+f)
+        RA_pnts.append(float(hdul[1].header['RA_PNT']))
+        DEC_pnts.append(float(hdul[1].header['DEC_PNT']))   
+        Tstart.append(Time(hdul[1].header['DATE-OBS']+'T'+hdul[1].header['TIME-OBS'], format='isot', scale='utc').mjd)
+        Tstop.append(Time(hdul[1].header['DATE-END']+'T'+hdul[1].header['TIME-END'], format='isot', scale='utc').mjd)
+        hdul.close()
+
+
+# In[5]:
+
+
+Coords_s=SkyCoord(RA,DEC,unit='degree')
+COORDS_pnts=SkyCoord(RA_pnts,DEC_pnts,unit='degree')
+seps=COORDS_pnts.separation(Coords_s).deg
+
+
+# In[6]:
+
+
+mask=np.where((seps<Radius) & (Tstart>T1) & (Tstop<T2))[0]
+OBSlist=[]
+for i in mask:
+    OBSlist.append(DL3_files[i])
+if(len(OBSlist)==0):
+    message='No data found'
+    raise RuntimeError('No data found')
+message
+
+
+# In[8]:
+
+
+cdec=cos(DEC*pi/180.)
+Npix=int(4*Radius/pixsize)+1
+RA_bins=np.linspace(RA-Radius/cdec,RA+Radius/cdec,Npix+1)
+DEC_bins=np.linspace(DEC-Radius,DEC+Radius,Npix+1)
+image=np.zeros((Npix,Npix))
+for f in OBSlist:
+    hdul=fits.open('data/'+f)
+    ev=hdul['EVENTS'].data
+    ev_ra=ev['RA']
+    ev_dec=ev['DEC']
+    ev_en=ev['ENERGY']
+    ev_time=ev['TIME']    
+    h=np.histogram2d(ev_ra,ev_dec,bins=[RA_bins,DEC_bins])
+    image+=h[0]
+    hdul.close()
+
+plt.imshow(np.flip(image,axis=1),extent=(RA_bins[-1],RA_bins[0],DEC_bins[0],DEC_bins[-1]),origin='lower')
+plt.colorbar()
+
+plt.xlabel('RA, degrees')
+plt.ylabel('DEC,degrees')
+plt.savefig('Image.png',format='png')
+
+
+# In[9]:
+
+
+# Create a new WCS object.  The number of axes must be set
+# from the start
+w = wcs.WCS(naxis=2)
+
+# Set up an "Airy's zenithal" projection
+# Vector properties may be set with Python lists, or Numpy arrays
+w.wcs.crpix = [Npix/2., Npix/2.]
+w.wcs.cdelt = np.array([pixsize/cdec, pixsize])
+w.wcs.crval = [RA, DEC]
+w.wcs.ctype = ["RA---AIR", "DEC--AIR"]
+w.wcs.set_pv([(2, 1, 45.0)])
+
+# Now, write out the WCS object as a FITS header
+header = w.to_header()
+
+# header is an astropy.io.fits.Header object.  We can use it to create a new
+# PrimaryHDU and write it to a file.
+hdu = fits.PrimaryHDU(image,header=header)
+hdu.writeto('Image.fits',overwrite=True)
+hdu=fits.open('Image.fits')
+im=hdu[0].data
+from astropy.wcs import WCS
+wcs = WCS(hdu[0].header)
+plt.subplot(projection=wcs)
+plt.imshow(im,  origin='lower')
+plt.grid(color='white', ls='solid')
+plt.xlabel('RA')
+plt.ylabel('Dec')
+
+
+# In[10]:
+
+
+bin_image = PictureProduct.from_file('Image.png')
+fits_image=ImageDataProduct.from_fits_file('Image.fits')
+
+
+# In[11]:
+
+
+picture = bin_image # http://odahub.io/ontology#ODAPictureProduct
+image = fits_image # http://odahub.io/ontology#Image
+
+
+# In[ ]:
+
+
+
+
+
+# In[ ]:
+
+
+_simple_outs, _oda_outs = [], []
+_galaxy_meta_data = {}
+_oda_outs.append(('out_Image_picture', 'picture_galaxy.output', picture))
+_oda_outs.append(('out_Image_image', 'image_galaxy.output', image))
+
+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'}
+
+for _outn, _outfn, _outv in _simple_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 _numpy_available and isinstance(_outv, np.ndarray):
+        with open(_galaxy_outfile_name, 'wb') as fd:
+            np.savez(fd, _outv)
+        _galaxy_meta_data[_outn] = {'ext': 'npz'}
+    else:
+        with open(_galaxy_outfile_name, 'w') as fd:
+            json.dump(_outv, fd)
+        _galaxy_meta_data[_outn] = {'ext': 'expression.json'}
+
+with open(os.path.join(_galaxy_wd, 'galaxy.json'), 'w') as fd:
+    json.dump(_galaxy_meta_data, fd)
+print("*** Job finished successfully ***")
+
diff --git a/tools/hess/Lightcurve.py b/tools/hess/Lightcurve.py
new file mode 100644
index 00000000..fdfe2951
--- /dev/null
+++ b/tools/hess/Lightcurve.py
@@ -0,0 +1,280 @@
+#!/usr/bin/env python
+# coding: utf-8
+
+# In[ ]:
+
+
+import json
+import os
+import shutil
+import sys
+
+try:
+    import numpy as np
+    _numpy_available = True
+except ImportError:
+    _numpy_available = False
+
+try:
+    from oda_api.json import CustomJSONEncoder
+except ImportError:
+    from json import JSONEncoder as CustomJSONEncoder
+
+_galaxy_wd = os.getcwd()
+
+
+# In[1]:
+
+
+import astropy.units as u
+from astropy.coordinates import SkyCoord
+
+import matplotlib.pyplot as plt
+from gammapy.data import DataStore
+from gammapy.makers import MapDatasetMaker
+from gammapy.makers.utils import make_theta_squared_table
+from gammapy.maps import Map, MapAxis, WcsGeom
+from astropy.io import fits
+import numpy as np
+from numpy import pi,cos,sin,sqrt,log10
+import os
+from astropy import wcs
+from astropy.io import fits
+from oda_api.data_products import PictureProduct
+from oda_api.data_products import ODAAstropyTable
+from oda_api.data_products import ImageDataProduct
+from astropy.time import Time
+
+
+# In[2]:
+
+
+if(os.path.exists('hess_dl3_dr1.tar.gz')==False):
+    get_ipython().system('wget https://zenodo.org/record/1421099/files/hess_dl3_dr1.tar.gz')
+    get_ipython().system('tar -zxvf hess_dl3_dr1.tar.gz')
+
+
+# In[8]:
+
+
+src_name='Crab' #http://odahub.io/ontology#AstrophysicalObject
+RA=83.628700  # http://odahub.io/ontology#PointOfInterestRA
+DEC = 22.014700 # http://odahub.io/ontology#PointOfInterestDEC
+T1='2003-10-09T13:16:00.0'# http://odahub.io/ontology#StartTime
+T2='2005-10-10T13:16:00.0' # http://odahub.io/ontology#EndTime
+Radius=2.5  #http://odahub.io/ontology#AngleDegrees
+R_s=0.2     #http://odahub.io/ontology#AngleDegrees
+Emin=100.    #http://odahub.io/ontology#Energy_GeV
+Emax=10000. #http://odahub.io/ontology#Energy_GeV
+NTbins=2 # http://odahub.io/ontology#Integer
+
+
+# In[ ]:
+
+
+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)
+
+
+# In[13]:
+
+
+Emin=Emin/1e3
+Emax=Emax/1e3
+T1=Time(T1, format='isot', scale='utc').mjd
+T2=Time(T2, format='isot', scale='utc').mjd
+message=''
+RA_pnts=[]
+DEC_pnts=[]
+DL3_files=[]
+OBSIDs=[]
+Tstart=[]
+Tstop=[]
+flist=os.listdir('data')
+for f in flist:
+    if(f[-7:]=='fits.gz'):
+        DL3_files.append(f)
+        OBSIDs.append(int(f[20:26]))
+        hdul=fits.open('data/'+f)
+        RA_pnts.append(float(hdul[1].header['RA_PNT']))
+        DEC_pnts.append(float(hdul[1].header['DEC_PNT']))   
+        Tstart.append(Time(hdul[1].header['DATE-OBS']+'T'+hdul[1].header['TIME-OBS'], format='isot', scale='utc').mjd)
+        Tstop.append(Time(hdul[1].header['DATE-END']+'T'+hdul[1].header['TIME-END'], format='isot', scale='utc').mjd)
+        hdul.close()
+
+
+# In[20]:
+
+
+Tbins=np.linspace(T1,T2,NTbins+1)
+Tmin=Tbins[:-1]
+Tmax=Tbins[1:]
+Tmean=(Tmin+Tmax)/2.
+Tbins
+
+
+# In[21]:
+
+
+Coords_s=SkyCoord(RA,DEC,unit='degree')
+COORDS_pnts=SkyCoord(RA_pnts,DEC_pnts,unit='degree')
+seps=COORDS_pnts.separation(Coords_s).deg
+
+
+# In[22]:
+
+
+mask=np.where((seps<Radius) & (Tstart>T1) & (Tstop<T2))[0]
+OBSlist=[]
+for i in mask:
+    OBSlist.append(DL3_files[i])
+if(len(OBSlist)==0):
+    message='No data found'
+    raise RuntimeError('No data found')
+message
+
+
+# In[23]:
+
+
+flux=np.zeros(NTbins)
+flux_err=np.zeros(NTbins)
+flux_b=np.zeros(NTbins)
+flux_b_err=np.zeros(NTbins)
+Expos=np.zeros(NTbins)
+for count,f in enumerate(OBSlist):
+    hdul=fits.open('data/'+f)
+    RA_pnt=hdul[1].header['RA_PNT']
+    DEC_pnt=hdul[1].header['DEC_PNT']
+    Texp=hdul[1].header['LIVETIME']
+    Trun_start=hdul[1].header['TSTART']
+    dRA=RA-RA_pnt
+    dDEC=DEC-DEC_pnt
+    RA_b=RA_pnt-dRA
+    DEC_b=DEC_pnt-dDEC
+    Coords_b=SkyCoord(RA_b,DEC_b,unit='degree')    
+    Coords_pnt=SkyCoord(RA_pnt,DEC_pnt,unit='degree')
+    dist=Coords_pnt.separation(Coords_s).deg
+
+    ev=hdul['EVENTS'].data
+    ev_ra=ev['RA']
+    ev_dec=ev['DEC']
+    ev_en=ev['ENERGY']
+    ev_time=(ev['TIME']-Trun_start)/86400.+Tbegs[count]
+    ev_coords=SkyCoord(ev_ra,ev_dec,unit='degree')
+    sep_s=ev_coords.separation(Coords_s).deg
+    sep_b=ev_coords.separation(Coords_b).deg
+
+    hdu=hdul['AEFF'].data
+    EEmin=hdu['ENERG_LO'][0]
+    EEmax=hdu['ENERG_HI'][0]
+    EE=sqrt(EEmin*EEmax)
+    EEbins=np.concatenate((EEmin,[EEmax[-1]]))
+    AA=hdu['EFFAREA'][0]+1e-10
+    Thmin=hdu['THETA_LO'][0]
+    Thmax=hdu['THETA_HI'][0]   
+    ind=np.argmin((Thmin-dist)**2)
+    AA=AA[ind]*Texp*1e4
+    mask=np.where((sep_s<R_s))
+    cts1=np.histogram2d(ev_time[mask],ev_en[mask],bins=[Tbins,EEbins])[0]
+    mask=(sep_b<R_s)
+    cts2=np.histogram2d(ev_time[mask],ev_en[mask],bins=[Tbins,EEbins])[0]
+    src=cts1-cts2
+    src_err=sqrt(cts1+cts2)
+    flux+=np.sum(src/AA,axis=1)
+    flux_err+=np.sum(src_err/AA,axis=1)
+    flux_b+=np.sum(cts2/AA,axis=1)
+    flux_b_err+=np.sum(sqrt(cts2)/AA,axis=1)
+    hdul.close()
+
+
+# In[24]:
+
+
+if(message==''):
+    plt.errorbar(Tmean,flux,yerr=flux_err,xerr=[Tmean-Tmin,Tmax-Tmean],linestyle='none',label='source')
+    plt.errorbar(Tmean,flux_b,yerr=flux_b_err,xerr=[Tmean-Tmin,Tmax-Tmean],linestyle='none',label='background')
+    plt.xlabel('Time, MJD')
+    plt.ylabel('Flux, cts/cm$^2$s')
+    plt.yscale('log')
+    ymin=min(min(flux-flux_err),min(flux_b-flux_b_err))
+    ymax=max(max(flux+flux_err),max(flux_b+flux_b_err))    
+    plt.ylim(ymin/2.,2*ymax)
+    plt.legend(loc='lower left')
+    plt.savefig('Lightcurve.png',format='png')
+
+
+# In[25]:
+
+
+if (message==''):
+    bin_image = PictureProduct.from_file('Lightcurve.png')
+    from astropy.table import Table
+    data=[Tmean,Tmin,Tmax,flux,flux_err,flux_b,flux_b_err]
+    names=('Tmean[MJD]','Tmin[MJD]','Tmax[MJD]','Flux[counts/cm2s]','Flux_error[counts/cm2s]','Background[counts/cm2s]','Background_error[counts/cm2s]')
+    lc = ODAAstropyTable(Table(data, names = names))
+
+
+# In[26]:
+
+
+picture = bin_image # http://odahub.io/ontology#ODAPictureProduct
+lightcurve_astropy_table = lc # http://odahub.io/ontology#ODAAstropyTable
+
+
+# In[ ]:
+
+
+
+
+
+# In[ ]:
+
+
+_simple_outs, _oda_outs = [], []
+_galaxy_meta_data = {}
+_oda_outs.append(('out_Lightcurve_picture', 'picture_galaxy.output', picture))
+_oda_outs.append(('out_Lightcurve_lightcurve_astropy_table', 'lightcurve_astropy_table_galaxy.output', lightcurve_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'}
+
+for _outn, _outfn, _outv in _simple_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 _numpy_available and isinstance(_outv, np.ndarray):
+        with open(_galaxy_outfile_name, 'wb') as fd:
+            np.savez(fd, _outv)
+        _galaxy_meta_data[_outn] = {'ext': 'npz'}
+    else:
+        with open(_galaxy_outfile_name, 'w') as fd:
+            json.dump(_outv, fd)
+        _galaxy_meta_data[_outn] = {'ext': 'expression.json'}
+
+with open(os.path.join(_galaxy_wd, 'galaxy.json'), 'w') as fd:
+    json.dump(_galaxy_meta_data, fd)
+print("*** Job finished successfully ***")
+
diff --git a/tools/hess/Spectrum.py b/tools/hess/Spectrum.py
new file mode 100644
index 00000000..6f899773
--- /dev/null
+++ b/tools/hess/Spectrum.py
@@ -0,0 +1,257 @@
+#!/usr/bin/env python
+# coding: utf-8
+
+# In[ ]:
+
+
+import json
+import os
+import shutil
+import sys
+
+try:
+    import numpy as np
+    _numpy_available = True
+except ImportError:
+    _numpy_available = False
+
+try:
+    from oda_api.json import CustomJSONEncoder
+except ImportError:
+    from json import JSONEncoder as CustomJSONEncoder
+
+_galaxy_wd = os.getcwd()
+
+
+# In[1]:
+
+
+import astropy.units as u
+from astropy.coordinates import SkyCoord
+import matplotlib.pyplot as plt
+from astropy.io import fits
+import numpy as np
+from numpy import pi,cos,sin,sqrt,log10
+import os
+
+from oda_api.data_products import PictureProduct
+from oda_api.data_products import ODAAstropyTable
+import os
+from astropy.time import Time
+if(os.path.exists('hess_dl3_dr1.tar.gz')==False):
+    get_ipython().system('wget https://zenodo.org/record/1421099/files/hess_dl3_dr1.tar.gz')
+    get_ipython().system('tar -zxvf hess_dl3_dr1.tar.gz')
+
+
+# In[2]:
+
+
+src_name='Crab' #http://odahub.io/ontology#AstrophysicalObject
+RA = 83.628700  # http://odahub.io/ontology#PointOfInterestRA
+DEC = 22.014700 # http://odahub.io/ontology#PointOfInterestDEC
+T1='2000-10-09T13:16:00.0'# http://odahub.io/ontology#StartTime
+T2='2022-10-10T13:16:00.0' # http://odahub.io/ontology#EndTime
+Radius=2.5  #http://odahub.io/ontology#AngleDegrees
+R_s=0.2     #http://odahub.io/ontology#AngleDegrees
+
+Emin=100.    #http://odahub.io/ontology#Energy_GeV
+Emax=10000. #http://odahub.io/ontology#Energy_GeV
+NEbins=20 # http://odahub.io/ontology#Integer
+
+
+# In[ ]:
+
+
+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)
+
+
+# In[3]:
+
+
+Emin=Emin/1e3
+Emax=Emax/1e3
+Ebins=np.logspace(log10(Emin),log10(Emax),NEbins+1)
+Ebins
+Emin=Ebins[:-1]
+Emax=Ebins[1:]
+Emean=sqrt(Emin*Emax)
+lgEmean=log10(Emean)
+
+
+# In[4]:
+
+
+T1=Time(T1, format='isot', scale='utc').mjd
+T2=Time(T2, format='isot', scale='utc').mjd
+message=''
+RA_pnts=[]
+DEC_pnts=[]
+DL3_files=[]
+OBSIDs=[]
+Tstart=[]
+Tstop=[]
+flist=os.listdir('data')
+for f in flist:
+    if(f[-7:]=='fits.gz'):
+        DL3_files.append(f)
+        OBSIDs.append(int(f[20:26]))
+        hdul=fits.open('data/'+f)
+        RA_pnts.append(float(hdul[1].header['RA_PNT']))
+        DEC_pnts.append(float(hdul[1].header['DEC_PNT']))   
+        Tstart.append(Time(hdul[1].header['DATE-OBS']+'T'+hdul[1].header['TIME-OBS'], format='isot', scale='utc').mjd)
+        Tstop.append(Time(hdul[1].header['DATE-END']+'T'+hdul[1].header['TIME-END'], format='isot', scale='utc').mjd)
+        hdul.close()
+
+
+# In[5]:
+
+
+Coords_s=SkyCoord(RA,DEC,unit='degree')
+COORDS_pnts=SkyCoord(RA_pnts,DEC_pnts,unit='degree')
+seps=COORDS_pnts.separation(Coords_s).deg
+
+
+# In[6]:
+
+
+mask=np.where((seps<Radius) & (Tstart>T1) & (Tstop<T2))[0]
+OBSlist=[]
+for i in mask:
+    OBSlist.append(DL3_files[i])
+if(len(OBSlist)==0):
+    message='No data found'
+    raise RuntimeError('No data found')
+message
+
+
+# In[9]:
+
+
+cts_s=np.zeros(NEbins)
+cts_b=np.zeros(NEbins)
+Expos=np.zeros(NEbins)
+for f in OBSlist:
+    hdul=fits.open('data/'+f)
+    RA_pnt=hdul[1].header['RA_PNT']
+    DEC_pnt=hdul[1].header['DEC_PNT']
+    Texp=hdul[1].header['LIVETIME']
+    dRA=RA-RA_pnt
+    dDEC=DEC-DEC_pnt
+    RA_b=RA_pnt-dRA
+    DEC_b=DEC_pnt-dDEC
+    Coords_b=SkyCoord(RA_b,DEC_b,unit='degree')    
+    Coords_pnt=SkyCoord(RA_pnt,DEC_pnt,unit='degree')
+    dist=Coords_pnt.separation(Coords_s).deg
+
+    ev=hdul['EVENTS'].data
+    ev_ra=ev['RA']
+    ev_dec=ev['DEC']
+    ev_en=ev['ENERGY']
+    ev_time=ev['TIME'] 
+    ev_coords=SkyCoord(ev_ra,ev_dec,unit='degree')
+    sep_s=ev_coords.separation(Coords_s).deg
+    sep_b=ev_coords.separation(Coords_b).deg
+
+    hdu=hdul['AEFF'].data
+    EEmin=hdu['ENERG_LO'][0]
+    EEmax=hdu['ENERG_HI'][0]
+    lgEE=log10(sqrt(EEmin*EEmax))
+    lgAA=log10(hdu['EFFAREA'][0]+1e-10)
+    Thmin=hdu['THETA_LO'][0]
+    Thmax=hdu['THETA_HI'][0]   
+    ind=np.argmin((Thmin-dist)**2)
+    Expos+=10**(np.interp(lgEmean,lgEE,lgAA[ind]))*Texp
+    mask=(sep_s<R_s)
+    cts_s+=np.histogram(ev_en[mask],bins=Ebins)[0]
+    mask=(sep_b<R_s)
+    cts_b+=np.histogram(ev_en[mask],bins=Ebins)[0]
+    hdul.close()
+
+
+# In[10]:
+
+
+flux=(cts_s-cts_b)/(Emax-Emin)*Emax*Emin/(Expos*1e4)
+flux_err=sqrt(cts_s+cts_b)/(Emax-Emin)*Emax*Emin/(Expos*1e4)
+plt.errorbar(Emean,flux,yerr=flux_err,xerr=[Emean-Emin,Emax-Emean])
+plt.xscale('log')
+plt.yscale('log')
+plt.xlabel('$E$, TeV')
+plt.ylabel('$E^2 dN/dE$, erg/(cm$^2$s)')
+plt.savefig('Spectrum.png',format='png')
+
+
+# In[11]:
+
+
+bin_image = PictureProduct.from_file('Spectrum.png')
+from astropy.table import Table
+data=[Emean,Emin,Emax,flux,flux_err,cts_s,cts_b,Expos*1e4]
+names=('Emean[TeV]','Emin[TeV]','Emax[TeV]','Flux[TeV/cm2s]','Flux_error[TeV/cm2s]','Cts_s','Cts_b','Exposure[cm2s]')
+spec = ODAAstropyTable(Table(data, names = names))
+
+
+# In[12]:
+
+
+picture_png = bin_image # http://odahub.io/ontology#ODAPictureProduct
+spectrum_astropy_table = spec # http://odahub.io/ontology#ODAAstropyTable
+
+
+# In[ ]:
+
+
+
+
+
+# In[ ]:
+
+
+_simple_outs, _oda_outs = [], []
+_galaxy_meta_data = {}
+_oda_outs.append(('out_Spectrum_picture_png', 'picture_png_galaxy.output', picture_png))
+_oda_outs.append(('out_Spectrum_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'}
+
+for _outn, _outfn, _outv in _simple_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 _numpy_available and isinstance(_outv, np.ndarray):
+        with open(_galaxy_outfile_name, 'wb') as fd:
+            np.savez(fd, _outv)
+        _galaxy_meta_data[_outn] = {'ext': 'npz'}
+    else:
+        with open(_galaxy_outfile_name, 'w') as fd:
+            json.dump(_outv, fd)
+        _galaxy_meta_data[_outn] = {'ext': 'expression.json'}
+
+with open(os.path.join(_galaxy_wd, 'galaxy.json'), 'w') as fd:
+    json.dump(_galaxy_meta_data, fd)
+print("*** Job finished successfully ***")
+
diff --git a/tools/hess/hess_astro_tool.xml b/tools/hess/hess_astro_tool.xml
new file mode 100644
index 00000000..e979c424
--- /dev/null
+++ b/tools/hess/hess_astro_tool.xml
@@ -0,0 +1,138 @@
+<tool id="hess_astro_tool" name="HESS" version="0.0.2+galaxy0" profile="23.0">
+  <requirements>
+    <requirement type="package" version="8.20.0">ipython</requirement>
+    <requirement type="package" version="5.3.4">astropy</requirement>
+    <requirement type="package" version="1.1">gammapy</requirement>
+    <requirement type="package" version="3.8.2">matplotlib</requirement>
+    <requirement type="package" version="0.4.6">astroquery</requirement>
+    <requirement type="package" version="1.11.4">scipy</requirement>
+    <requirement type="package" version="1.2.5">oda-api</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.14.2">nbconvert</requirement>
+    <requirement type="package" version="1.20.3">wget</requirement>
+  </requirements>
+  <command detect_errors="exit_code">ipython '$__tool_directory__/${_data_product._selector}.py'</command>
+  <configfiles>
+    <inputs name="inputs" filename="inputs.json" />
+  </configfiles>
+  <inputs>
+    <conditional name="_data_product">
+      <param name="_selector" type="select" label="Data Product">
+        <option value="Image" selected="true">Image</option>
+        <option value="Spectrum" selected="false">Spectrum</option>
+        <option value="Lightcurve" selected="false">Lightcurve</option>
+      </param>
+      <when value="Image">
+        <param name="src_name" type="text" value="Crab" label="src_name" />
+        <param name="RA" type="float" value="83.6287" label="RA Units: deg" />
+        <param name="DEC" type="float" value="22.0147" label="DEC Units: deg" />
+        <param name="T1" type="text" value="2000-10-09T13:16:00.0" label="T1" />
+        <param name="T2" type="text" value="2022-10-10T13:16:00.0" label="T2" />
+        <param name="Radius" type="float" value="2.5" label="Radius Units: deg" />
+        <param name="pixsize" type="float" value="0.1" label="pixsize Units: deg" />
+        <param name="Emin" type="float" value="100.0" label="Emin Units: GeV" />
+        <param name="Emax" type="float" value="10000.0" label="Emax Units: GeV" />
+      </when>
+      <when value="Spectrum">
+        <param name="src_name" type="text" value="Crab" label="src_name" />
+        <param name="RA" type="float" value="83.6287" label="RA Units: deg" />
+        <param name="DEC" type="float" value="22.0147" label="DEC Units: deg" />
+        <param name="T1" type="text" value="2000-10-09T13:16:00.0" label="T1" />
+        <param name="T2" type="text" value="2022-10-10T13:16:00.0" label="T2" />
+        <param name="Radius" type="float" value="2.5" label="Radius Units: deg" />
+        <param name="R_s" type="float" value="0.2" label="R_s Units: deg" />
+        <param name="Emin" type="float" value="100.0" label="Emin Units: GeV" />
+        <param name="Emax" type="float" value="10000.0" label="Emax Units: GeV" />
+        <param name="NEbins" type="integer" value="20" label="NEbins" />
+      </when>
+      <when value="Lightcurve">
+        <param name="src_name" type="text" value="Crab" label="src_name" />
+        <param name="RA" type="float" value="83.6287" label="RA Units: deg" />
+        <param name="DEC" type="float" value="22.0147" label="DEC Units: deg" />
+        <param name="T1" type="text" value="2003-10-09T13:16:00.0" label="T1" />
+        <param name="T2" type="text" value="2005-10-10T13:16:00.0" label="T2" />
+        <param name="Radius" type="float" value="2.5" label="Radius Units: deg" />
+        <param name="R_s" type="float" value="0.2" label="R_s Units: deg" />
+        <param name="Emin" type="float" value="100.0" label="Emin Units: GeV" />
+        <param name="Emax" type="float" value="10000.0" label="Emax Units: GeV" />
+        <param name="NTbins" type="integer" value="2" label="NTbins" />
+      </when>
+    </conditional>
+  </inputs>
+  <outputs>
+    <data label="${tool.name} -&gt; Image picture" name="out_Image_picture" format="auto" from_work_dir="picture_galaxy.output">
+      <filter>_data_product['_selector'] == 'Image'</filter>
+    </data>
+    <data label="${tool.name} -&gt; Image image" name="out_Image_image" format="auto" from_work_dir="image_galaxy.output">
+      <filter>_data_product['_selector'] == 'Image'</filter>
+    </data>
+    <data label="${tool.name} -&gt; Spectrum picture_png" name="out_Spectrum_picture_png" format="auto" from_work_dir="picture_png_galaxy.output">
+      <filter>_data_product['_selector'] == 'Spectrum'</filter>
+    </data>
+    <data label="${tool.name} -&gt; Spectrum spectrum_astropy_table" name="out_Spectrum_spectrum_astropy_table" format="auto" from_work_dir="spectrum_astropy_table_galaxy.output">
+      <filter>_data_product['_selector'] == 'Spectrum'</filter>
+    </data>
+    <data label="${tool.name} -&gt; Lightcurve picture" name="out_Lightcurve_picture" format="auto" from_work_dir="picture_galaxy.output">
+      <filter>_data_product['_selector'] == 'Lightcurve'</filter>
+    </data>
+    <data label="${tool.name} -&gt; Lightcurve lightcurve_astropy_table" name="out_Lightcurve_lightcurve_astropy_table" format="auto" from_work_dir="lightcurve_astropy_table_galaxy.output">
+      <filter>_data_product['_selector'] == 'Lightcurve'</filter>
+    </data>
+  </outputs>
+  <tests>
+    <test expect_num_outputs="2">
+      <conditional name="_data_product">
+        <param name="_selector" value="Image" />
+        <param name="src_name" value="Crab" />
+        <param name="RA" value="83.6287" />
+        <param name="DEC" value="22.0147" />
+        <param name="T1" value="2000-10-09T13:16:00.0" />
+        <param name="T2" value="2022-10-10T13:16:00.0" />
+        <param name="Radius" value="2.5" />
+        <param name="pixsize" value="0.1" />
+        <param name="Emin" value="100.0" />
+        <param name="Emax" value="10000.0" />
+      </conditional>
+      <assert_stdout>
+        <has_text text="*** Job finished successfully ***" />
+      </assert_stdout>
+    </test>
+    <test expect_num_outputs="2">
+      <conditional name="_data_product">
+        <param name="_selector" value="Spectrum" />
+        <param name="src_name" value="Crab" />
+        <param name="RA" value="83.6287" />
+        <param name="DEC" value="22.0147" />
+        <param name="T1" value="2000-10-09T13:16:00.0" />
+        <param name="T2" value="2022-10-10T13:16:00.0" />
+        <param name="Radius" value="2.5" />
+        <param name="R_s" value="0.2" />
+        <param name="Emin" value="100.0" />
+        <param name="Emax" value="10000.0" />
+        <param name="NEbins" value="20" />
+      </conditional>
+      <assert_stdout>
+        <has_text text="*** Job finished successfully ***" />
+      </assert_stdout>
+    </test>
+    <test expect_num_outputs="2">
+      <conditional name="_data_product">
+        <param name="_selector" value="Lightcurve" />
+        <param name="src_name" value="Crab" />
+        <param name="RA" value="83.6287" />
+        <param name="DEC" value="22.0147" />
+        <param name="T1" value="2003-10-09T13:16:00.0" />
+        <param name="T2" value="2005-10-10T13:16:00.0" />
+        <param name="Radius" value="2.5" />
+        <param name="R_s" value="0.2" />
+        <param name="Emin" value="100.0" />
+        <param name="Emax" value="10000.0" />
+        <param name="NTbins" value="2" />
+      </conditional>
+      <assert_stdout>
+        <has_text text="*** Job finished successfully ***" />
+      </assert_stdout>
+    </test>
+  </tests>
+</tool>
\ No newline at end of file