Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

New Sources added #242

Open
wants to merge 16 commits into
base: master
Choose a base branch
from
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Binary file added .DS_Store
Binary file not shown.
54 changes: 54 additions & 0 deletions HESS_J1825/Index_vs_flux.txt
Original file line number Diff line number Diff line change
@@ -0,0 +1,54 @@
# %ECSV 0.9
# ---
# datatype:
# - {name: boxid, datatype: string, format: '%.2f'}
# - {name: dist_deg, unit: deg, datatype: float32, format: '%.2f'}
# - {name: dist_pc, unit: pc, datatype: float32, format: '%.2f'}
# - {name: gamma, datatype: float32, format: '%.2f'}
# - {name: gamma_stat, datatype: float32, format: '%.2f'}
# - {name: gamma_sys, datatype: float32, format: '%.2f'}
# - {name: dnde_1_5, unit: 1 / (cm2 TeV s), datatype: float32, format: '%2.4g'}
# - {name: dnde_1_5_stat_err, unit: 1 / (cm2 TeV s), datatype: float32, format: '%2.4g'}
# - {name: dnde_1_5_sys_err, unit: 1 / (cm2 TeV s), datatype: float32, format: '%2.4g'}
# - {name: sigma, datatype: float32, format: '%2.4g'}
# meta: !!omap
# schema: astropy-2.0
#boxid dist_deg dist_pc gamma gamma_stat gamma_sys dnde_1_5 dnde_1_5_stat_err dnde_1_5_sys_err sigma
1 0.74 51 2.32 0.05 0.03 4.3e-13 0.2e-13 0.2e-13 28.8
2 0.58 41 2.21 0.04 0.03 5.9e-13 0.2e-13 0.3e-13 34.2
3 0.52 36 2.24 0.06 0.03 4.3e-13 0.2e-13 0.2e-13 24.6
4 0.58 41 2.29 0.04 0.03 6.0e-13 0.2e-13 0.2e-13 38.8
5 0.36 26 2.13 0.03 0.03 9.6e-13 0.3e-13 0.4e-13 51.4
6 0.26 18 2.09 0.04 0.02 6.9e-13 0.3e-13 0.3e-13 35.05
7 0.52 36 2.20 0.04 0.02 5.4e-13 0.2e-13 0.2e-13 33.7
8 0.26 18 2.06 0.03 0.01 10.4e-13 0.3e-13 0.3e-13 53.6
9 0. 0 2.00 0.03 0.02 8.7e-13 0.3e-13 0.3e-13 41.6
10 0.58 41 2.29 0.07 0.05 2.9e-13 0.2e-13 0.2e-13 20.7
11 0.37 26 2.20 0.05 0.03 4.3e-13 0.2e-13 0.1e-13 26.4
12 0.26 18 2.06 0.06 0.04 4.0e-13 0.2e-13 0.2e-13 20.5
a 1.52 106 3.3 0.2 0.2 0.5e-13 0.1e-13 0.01e-13 8.8
b 1.40 98 3.2 0.2 0.8 0.6e-13 0.2e-13 0.4e-13 6.9
c 1.33 93 3.3 0.3 0.3 0.5e-13 0.2e-13 0.3e-13 6.9
d 1.47 103 3.1 0.3 0.4 0.5e-13 0.2e-13 0.01e-13 6.2
e 1.3 91 2.90 0.13 0.05 1.2e-13 0.2e-13 0.1e-13 12.8
f 1.16 81 2.8 0.1 0.1 1.3e-13 0.2e-13 0.2e-13 13.3
g 1.07 75 2.7 0.1 0.1 1.3e-13 0.2e-13 0.2e-13 11.1
h 1.04 73 2.9 0.2 0.1 0.8e-13 0.2e-13 0.1e-13 8.2
i 1.3 91 3.1 0.2 0.1 0.7e-13 0.2e-13 0.2e-13 9.4
j 1.10 77 2.61 0.09 0.07 2.1e-13 0.2e-13 0.2e-13 17.7
k 0.94 65 2.43 0.07 0.04 2.8e-13 0.2e-13 0.1e-13 20.4
l 0.82 57 2.40 0.07 0.03 3.1e-13 0.2e-13 0.2e-13 21.0
m 0.78 56 2.40 0.10 0.05 2.2e-13 0.2e-13 0.1e-13 14.5
n 0.82 57 2.45 0.15 0.3 1.3e-13 0.2e-13 0.1e-13 8.7
o 1.40 98 2.8 0.2 0.1 0.9e-13 0.2e-13 0.1e-13 9.7
p 1.16 81 2.72 0.11 0.09 1.6e-13 0.2e-13 0.2e-13 15.1
q 0.94 65 2.33 0.07 0.02 3.0e-13 0.2e-13 0.2e-13 21.2
r 0.58 41 2.59 0.1 0.03 1.9e-13 0.2e-13 0.2e-13 13.8
s 1.07 75 2.8 0.1 0.1 1.1e-13 0.2e-13 0.2e-13 13.7
t 0.82 57 2.3 0.07 0.02 3.2e-13 0.2e-13 0.2e-13 22.4
u 0.78 56 2.28 0.08 0.07 2.6e-13 0.2e-13 0.3e-13 19.2
v 0.82 57 2.47 0.13 0.06 1.4e-13 0.2e-13 0.1e-13 12.4
w 0.94 65 2.45 0.16 0.1 1.3e-13 0.2e-13 0.2e-13 10.7
x 0.74 51 2.71 0.16 0.2 1.0e-13 0.2e-13 0.2e-13 10.6
y 0.58 41 2.5 0.2 0.1 1.2e-13 0.2e-13 0.1e-13 9.8
z 0.52 36 2.6 0.2 0.09 1.0e-13 0.2e-13 0.1e-13 8.1
154 changes: 154 additions & 0 deletions HESS_J1825/MAP_CREATOR.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,154 @@
from gammapy.maps import Map, MapAxis
from astropy.coordinates import SkyCoord
#from gammapy.spectrum.models import PowerLaw,PowerLaw2,ExponentialCutoffPowerLaw
import numpy as np
from astropy import units as u
from astropy.io import fits
import matplotlib.pyplot as plt
from scipy import interpolate
from gammapy.maps import Map, WcsGeom, WcsNDMap
from gammapy.modeling.models import DiskSpatialModel, PointSpatialModel, ExpCutoffPowerLawSpectralModel
import gammalib

#######
def pow(x,N,indice,eref):
return N*(x/eref)**(-1.*indice)

def expcutoff(x,N,indice,ecut,eref):
return N*(x/eref)**(-1.*indice) * np.exp(-1.0*x/ecut)
########

### read the file with the spatial extent
rad = np.genfromtxt('radial_extent.txt')
energy = rad[:,0]*1e6
dist = rad[:,1] * 1.357 #the factor 1.357 relates the radius at 1/e with the radius at 90% (0.9/(1/e)**1/3)

rad_interpolate = interpolate.interp1d(energy, dist, fill_value="extrapolate")

### read the file with the spectral parameters and estimate the normalization
spec = np.genfromtxt('Index_vs_flux.txt')
index = spec[:,3]
e_index = spec[:,4]
ang_dist = spec[:,1]
flux_1_5TeV = spec[:,6]
e_flux_1_5TeV = spec[:,7]

### I've estimated a function to model the index and fluxes dependecies, based on Fig. 9 of HESS paper
def fluxes(x):
a = 1.2 * 1.08686e-12
b = 0.0326925
c = 0.528341
return a * (x + 1)**b * np.exp(-1.0*x/c)

def indexes(x):
a = 0.697262
b = 1.87344
return a*x+b

# Define energy ranges and a cut-off
emin=1e6 #TeV
emax=5e6 #TeV
E_ref = 1e6 #MeV
Ecut = 19e6 # cutoff in MeV

### interpolate the function of index, fluxes and normalizations
d = np.linspace(0.0,1.55,100)

### Position of the pulsar
ra_0_psr = 276.554417
dec_0_psr = -13.580028
position_psr = SkyCoord(ra_0_psr, dec_0_psr, frame='icrs', unit='deg')

### Position of the map center
position = SkyCoord(276.4, -13.966667, frame='icrs', unit='deg')

### Energy axes and size
energy_axis = MapAxis.from_bounds(0.2e5, 250e6, 40, interp='log', name='energy', unit='MeV')

enlarger = 1.2
r_0 = rad_interpolate(energy_axis.center) * u.deg * enlarger
r_0 = np.where(r_0>=0, r_0, 0) * u.deg

### Fill the index and normalization at the positions of the pixels
binsize=0.02
eccentricity = 0.69
ang = 45.
width = 3.

m_wcs = WcsGeom.create(binsz=binsize, coordsys="CEL",
skydir=position_psr,
width=width*u.deg,
axes=[energy_axis])

model = DiskSpatialModel(lon_0=ra_0_psr*u.deg, lat_0=dec_0_psr*u.deg, r_0="1.1 deg", e=eccentricity, phi=ang * u.deg, frame="icrs")

coords = m_wcs.get_coord()
coord_map = SkyCoord(coords.lon, coords.lat, frame='icrs')
separation = position_psr.separation(coord_map)

new_indexes = indexes(separation.value)

N0 = (fluxes(separation.value) * (1-new_indexes) /
((emax/E_ref)**(1-new_indexes) - (emin/E_ref)**(1-new_indexes))) / (0.26 * np.pi / 180.)**2

### Create the map
skymap = Map.create(width=width*u.deg, axes=[energy_axis],
binsz=binsize, coordsys="CEL",
skydir=position_psr)

for idx, energy in enumerate(energy_axis.center):

if r_0[idx] > 0:
model = DiskSpatialModel(lon_0=ra_0_psr*u.deg, lat_0=dec_0_psr*u.deg,
r_0=r_0[idx],
e=eccentricity,
phi=ang * u.deg,
frame="icrs")
else:
model = PointSpatialModel(lon_0=ra_0_psr*u.deg, lat_0=dec_0_psr*u.deg,
frame='icrs')

skymap.data[idx] = model.evaluate_geom(skymap.geom.to_image())

skymap.data = np.where(skymap.data>0,1.,0)

### Fill the map with the flux values
for j in np.arange(len(energy_axis.edges)-1):
flx = (expcutoff((energy_axis.edges.value[j]/1e6),N0[j],new_indexes[j],Ecut/1e6,E_ref/1e6)
* (energy_axis.edges.value[j+1]-energy_axis.edges.value[j])/1e6)

flx = flx / (energy_axis.edges.value[j+1]-energy_axis.edges.value[j])
skymap.data[j] *= flx

### Cut the ellipse in order to get an asymettric shape
for j in np.arange(len(energy_axis.edges)-1):
for idx1 in np.arange(0,int(width/binsize)):
for idx2 in (np.arange(0,int(width/binsize))[::-1]):
if (idx2+1.0)/(idx1+1.0) < (np.tan((125-ang)*np.pi/180.) - 7/binsize/(idx1+1.0)): #+ 80/idx1+1.0:
skymap.data[j,idx1,idx2] *= 0
else:
continue

### Smooth the map
skymap = skymap.smooth(width=0.08 * u.deg, kernel="gauss")

### Save the map
skymap.write('hessj1825_cube.fits', hdu='Primary', overwrite='True')

hdul = fits.open('hessj1825_cube.fits')
hdul[0].header.set('BUNIT', 'photon/cm2/s/MeV/sr', 'Photon flux')
hdul[0].verify('fix')
hdul[1].name='ENERGIES'
hdul[1].header['TTYPE2']='Energy'
hdul[1].verify('fix')
hdul.writeto('hessj1825_map.fits',overwrite=True)

### create gammalib model
spatial = gammalib.GModelSpatialDiffuseCube(gammalib.GFilename('hessj1825_map.fits'))
spectral = gammalib.GModelSpectralConst(1.)
model = gammalib.GModelSky(spatial,spectral)
## fill to model container and write to disk
models = gammalib.GModels()
models.append(model)
models.save('hessj1825.xml')

11 changes: 11 additions & 0 deletions HESS_J1825/hessj1825.xml
Original file line number Diff line number Diff line change
@@ -0,0 +1,11 @@
<?xml version="1.0" encoding="UTF-8" standalone="no"?>
<source_library title="source library">
<source name="" type="DiffuseSource">
<spectrum type="Constant">
<parameter name="Normalization" value="1" error="0" scale="1" min="0" max="1000" free="1" />
</spectrum>
<spatialModel type="DiffuseMapCube" file="hessj1825_map.fits">
<parameter name="Normalization" value="1" scale="1" min="0.001" max="1000" free="0" />
</spatialModel>
</source>
</source_library>
Binary file added HESS_J1825/hessj1825_map.fits
Binary file not shown.
9 changes: 9 additions & 0 deletions HESS_J1825/radial_extent.txt
Original file line number Diff line number Diff line change
@@ -0,0 +1,9 @@
0.15176073466832585 0.6770880692315371
0.406158598837698 0.778267746075181
0.8463231720621138 0.7360977391991531
1.6774040264261578 0.6524032698070652
3.4952501326754755 0.47579443140094124
6.700187503509583 0.3915091671561024
12.84386263589743 0.27005462091056004
25.88472040318299 0.19694948009891955
47.99095809323868 0.1449740670372633
31 changes: 31 additions & 0 deletions gammapy-0.10-environment.yml
Original file line number Diff line number Diff line change
@@ -0,0 +1,31 @@
# Conda environment for Gammapy version 0.10
# Install: conda env create -f gammapy-0.10-environment.yml
# Activate: conda activate gammapy-0.10
# Deactivate: conda deactivate

name: gammapy-0.10

channels:
- conda-forge
- sherpa

dependencies:
- gammapy=0.10
- python=3.6
- ipython=7.2
- jupyter=1.0
- jupyterlab=0.35
- pyyaml=3.13
- click=7.0
- numpy=1.16
- scipy=1.2
- pandas=0.24
- matplotlib>=2,<3
- astropy=3.1
- regions=0.3
- reproject=0.4
- uncertainties=3.0
- iminuit=1.3
- sherpa=4.10
- healpy=1.11
- naima=0.8
Binary file added input/.DS_Store
Binary file not shown.
Binary file added input/data/2009/.DS_Store
Binary file not shown.
Binary file added input/data/2009/2009ApJ...703L...6A/.DS_Store
Binary file not shown.
2 changes: 1 addition & 1 deletion input/data/2009/2009ApJ...703L...6A/info.yaml
Original file line number Diff line number Diff line change
Expand Up @@ -7,5 +7,5 @@ datasets:
- tev-000149-sed.ecsv

data_entry:
status: complete
status: uncomplete
reviewed: no
Binary file added input/data/2011/.DS_Store
Binary file not shown.
10 changes: 10 additions & 0 deletions input/data/2011/2011arXiv1110.6368B/info.yaml
Original file line number Diff line number Diff line change
@@ -0,0 +1,10 @@
reference_id: 2011arXiv1110.6368B
source_id: [139]
telescope: [magic]

datasets:
- tev-000139.yaml

data_entry:
status: uncomplete
reviewed: no
15 changes: 15 additions & 0 deletions input/data/2011/2011arXiv1110.6368B/tev-000139.yaml
Original file line number Diff line number Diff line change
@@ -0,0 +1,15 @@
source_id: 139
reference_id: 2011arXiv1110.6368B
telescope: magic

data:
livetime: 1.5h
significance: 7.6

pos:
ra: {val: 20h01m13.5s}
dec: {val: 43d53m02.8s}

morph:
type: point

Binary file added input/data/2015/.DS_Store
Binary file not shown.
Binary file added input/data/2015/2015MNRAS.446..217A/.DS_Store
Binary file not shown.
2 changes: 1 addition & 1 deletion input/data/2015/2015MNRAS.446..217A/info.yaml
Original file line number Diff line number Diff line change
Expand Up @@ -6,5 +6,5 @@ datasets:
- tev-000005.yaml

data_entry:
status: missing
status: uncomplete
reviewed: no
22 changes: 20 additions & 2 deletions input/data/2015/2015MNRAS.446..217A/tev-000005.yaml
Original file line number Diff line number Diff line change
@@ -1,5 +1,23 @@
source_id: 5
reference_id: 2015MNRAS.446..217A
reference_id: 2015MNRAS.446..217A
telescope: magic

#TODO: Add data
data:
livetime: 19.7h
significance: 5.5

pos:
ra: {val: 00h35m16.8s, err: 0h0m7.2s}
dec: {val: +59d47m24.0s, err: 0h1m12.0s}

morph:
type: point

spec:
model:
type: pl
parameters:
norm: {val: 2.0, err: 0.50, err_sys: 0.50, scale: 1e-11, unit: cm-2 s-1 TeV-1}
index: {val: 3.8, err: 0.7, err_sys: 0.3}
e_ref: {val: 0.25, unit: TeV}
erange: {min: 0.125, max: 0.5, unit: TeV}
Binary file added input/data/2018/2018A&A...612A...1H/.DS_Store
Binary file not shown.
18 changes: 18 additions & 0 deletions input/data/2018/2018A&A...612A...1H/info.yaml
Original file line number Diff line number Diff line change
@@ -0,0 +1,18 @@
reference_id: 2018A&A...612A...1H
source_id: [50,122,127,129,160,162,163,164,170]
telescope: [hess]

datasets:
- tev-000050.yaml
- tev-000122.yaml
- tev-000127.yaml
- tev-000129.yaml
- tev-000160.yaml
- tev-000162.yaml
- tev-000163.yaml
- tev-000164.yaml
- tev-000170.yaml

data_entry:
status: uncomplete
reviewed: no
24 changes: 24 additions & 0 deletions input/data/2018/2018A&A...612A...1H/tev-000050.yaml
Original file line number Diff line number Diff line change
@@ -0,0 +1,24 @@
source_id: 50
reference_id: 2018A&A...612A...1H
telescope: hess

data:
livetime: 51.7h
significance: 10.2

pos:
ra: {val: 11h19m04.03s}
dec: {val: -61d27m07.8s}

morph:
type: gauss
sigma: {val: 5.9m, err: 0.9m}

spec:
model:
type: pl
parameters:
norm: {val: 7.96, err: 0.75, err_sys: 0.0, scale: 1e-13, unit: cm-2 s-1 TeV-1}
index: {val: 2.64, err: 0.12, err_sys: 0.2}
e_ref: {val: 1.27, unit: TeV}
erange: {min: 0.42, max: 68.13, unit: TeV}
Loading