In [1]:
import yaml
import matplotlib.pyplot
import numpy as np
from astropy.table import Table, join
from haloanalysis.sed import SED
from haloanalysis.utils import load_source_rows
from ebltable.tau_from_model import OptDepth
from fermipy.spectrum import PLExpCutoff
from collections import OrderedDict
from haloanalysis import model
import matplotlib.patches as mpatches
from os import path
In [2]:
%matplotlib inline
In [29]:
cfile = '../config/igmf_analysis/test.yaml'
cfile = '../config/igmf_analysis/elmag_1e1.yaml'
#cfile = '../config/igmf_analysis/analytic_1e7.yaml'
#cfile = '../config/igmf_analysis/elmag_1e4.yaml'
cfile = '../config/igmf_analysis/elmag_1e7.yaml'
In [30]:
config = yaml.load(open(cfile))
In [31]:
outfile = path.join(config['outdir'],'fit_igmf_th_jet{0[th_jet]:.0f}_tmax{0[tmax]:.0e}_lp_{0[kind]:s}_combined.fits'.format(config))
In [32]:
if config['int_spec'] == 'PLExpCutOff':
fint = sp.PLExpCutoff # intrinsic spectrum
elif config['int_spec'] == 'LogParabolaExpCutoff':
fint = model.LogParabolaExpCutoff # intrinsic spectrum
scale = 1e3
In [33]:
results = Table.read(outfile)
tab_sed_tev = Table.read(config['cat_tev'])
tables = [Table.read(config['cat'],'CATALOG'),
Table.read(config['cat'],'LIKELIHOOD'),
Table.read(config['cat'],'SED')]
tab_casc = join(tables[0],tables[1])
tab_casc = join(tab_casc,tables[2])
tab_pars = Table.read(config['cat'],hdu='SCAN_PARS')
tab_ebounds = Table.read(config['cat'],hdu='EBOUNDS')
In [34]:
results
Out[34]:
print the best fit results
In [35]:
idL, idB = 0,0
for isource, fit_pars in enumerate(results['src_fit_pars']):
print results['SOURCE_FULL'][isource], fit_pars[idL,idB]
In [36]:
tau = OptDepth.readmodel(model = config['ebl_model'])
EMeV = np.logspace(3,8,200)
In [37]:
cp = plt.cm.afmhot
idL = 8 # coh length = 1 Mpc
In [38]:
for isource in range(results['SOURCE_FULL'].shape[0]):
ax = plt.subplot(111)
# select the source
name = results['SOURCE_FULL'][isource]
rows_sed_tev = load_source_rows(tab_sed_tev, [name], key='SOURCE_FULL')
cat_names = [ row['3FGL_NAME'] for row in rows_sed_tev ]
cat_names = np.unique(np.array(cat_names))
rows_sed_gev = load_source_rows(tab_casc, cat_names, key='name_3fgl')
rows_casc = load_source_rows(tab_casc, cat_names, key='name_3fgl')
# load the source SEDs
sed_Fermi = SED.create_from_row2(rows_sed_gev, tab_ebounds)
sed_IACT = SED.create_from_row(rows_sed_tev)
efct_Fermi = sed_Fermi.ectr**2/sed_Fermi.ewidth*1.602177e-06
efct_IACT = sed_IACT.ectr**2/sed_IACT.ewidth*1.602177e-06
# plot the SED
sed_Fermi.plot(label = 'Fermi-LAT PS',
color = plt.cm.Blues_r(0.2),
mec = plt.cm.Blues_r(0.2),
zorder = 5)
sed_IACT.plot(label = 'IACT PS',
color = plt.cm.Blues_r(0.5),
mec = plt.cm.Blues_r(0.5),
marker = 's', zorder = 5)
for j,i in enumerate(([0,8,10,12,14,16])): # B field models
icol = j/(6. + 4.) + 0.2
if not isource: print icol
plt.plot(sed_Fermi.ectr/1E6,
efct_Fermi*results['casc_flux'][isource][idL,i],
label='$B = 10^{{{0:.1f}}}$ G'.format(
results['igmf'][isource,0,i]),
color = cp(icol),
lw = 1.5, ls = '-',
)
# intrinsic spectrum
fplot = fint(results['src_fit_pars'][isource,idL,i], scale = scale)
plt.plot(EMeV / 1e6,
fplot.e2dfde(EMeV)*1.602177e-06 * \
np.exp(-1. * tau.opt_depth(results['redshift'][isource], EMeV / 1e6)),
color = cp(icol),
ls = '--', lw = 1)
# Halo Flux UL
if not i:
plt.errorbar(sed_Fermi.ectr/1E6,
results['halo_flux_ul'][isource,idL,i]*efct_Fermi,
ls = 'none',
marker = 'v',
color = '0.7',
mec = '0.7',
uplims = True, zorder = -1, lw = 2,
alpha = 1.,
label = 'Halo UL for $B = 10^{{{0:.1f}}}$ G'.format(
results['igmf'][isource,0,i]) if not i else ''
)
title = '{0:s}; '.format(name.split('_')[0])
title += '$\lambda = 10^{{{0:.0f}}}\,$Mpc'.format(results['lcoh'][isource,idL,0])
title += ', $t_\mathrm{{max}} = 10^{0:.0f}$ yrs, $\\theta_\mathrm{{jet}} = {1:.1f}^\circ$'.format(
np.log10(config['tmax']), config['th_jet'])
#title += '\n{0:s}'.format(" ".join(name.split('_')))
fname = 'sed_{0:s}_idL{1:n}_tmax{2[tmax]:.0e}_thj{2[th_jet]:.0f}_{2[kind]:s}'.format(name, idL, config)
plt.legend(loc = 'lower center', fontsize = 'x-small', ncol = 3,
title = title, fancybox = True, framealpha = 0.9)
plt.gca().get_legend().get_title().set_fontsize('x-small')
plt.gca().set_xscale('log')
plt.gca().set_yscale('log')
plt.gca().set_ylim(1e-14,3e-11)
plt.gca().set_xlim(1e-3,1e1)
plt.xlabel('Energy (TeV)')
plt.ylabel('$E^2 dN/dE\,(\mathrm{erg}\,\mathrm{s}^{-1}\mathrm{cm}^{-2})$')
plt.savefig('./plots/' + fname + '.png', format = 'png', dpi = 200)
plt.savefig('./plots/' + fname + '.pdf', format = 'pdf')
plt.close()
print isource, fname
In [13]:
isource = 0
In [14]:
from myplot.funcs import *
In [15]:
cmap = plt.cm.afmhot
Get the bin boundaries
In [16]:
igmf = results['igmf'][isource]*np.ones(results['dloglike'][isource].shape)
lcoh = results['lcoh'][isource]*np.ones(results['dloglike'][isource].shape)
dl = results['lcoh'][isource][1] - results['lcoh'][isource][0]
lcoh_edges = np.linspace(results['lcoh'][isource][0] - dl[0] / 2.,
results['lcoh'][isource][-1] + dl[0] / 2.,
results['lcoh'][isource].shape[0] + 1)
dB = results['igmf'][isource][0,1] - results['igmf'][isource][0,0]
igmf_edges = np.linspace(results['igmf'][isource][0,0] - dB / 2.,
results['igmf'][isource][0,-1] + dB / 2.,
results['igmf'][isource].shape[1] + 1)
bb,ll = np.meshgrid(igmf_edges,lcoh_edges)
print bb.shape, ll.shape
In [17]:
cp = plt.cm.inferno
Get min and max likelihoods for color scale, check which spectra don't give a limit and exclude sources that give limits that are more than 10% different from each other (from excluded area)
In [18]:
no_lim = []
max_diff = 100. # max diff allowed between limits from different spectra
exclAreas = OrderedDict()
exclAreas_fhes = OrderedDict()
combined_dl_src = OrderedDict()
for isource, name in enumerate(results['SOURCE_FULL']):
dl = 2. * np.min(results['dloglike'][isource])
exclArea = np.sum(2. * results['dloglike'][isource] < -6.18) \
/ float(np.prod(results['dloglike'][isource].shape))
if not exclArea:
print 'No exclusion from ', name, results['name'][isource]
no_lim.append(results['name'][isource])
else:
# exclude 1426 by hand
if name.find('1426') >= 0:
no_lim.append(results['name'][isource])
continue
if not results['name'][isource] in exclAreas_fhes:
exclAreas_fhes[results['name'][isource]] = exclArea
elif np.abs(1. - exclAreas_fhes[results['name'][isource]] / exclArea) > max_diff:
print 'Large difference in excluded area for different spectra, removing ', name
no_lim.append(results['name'][isource])
continue
exclAreas[name] = exclArea
if not results['name'][isource] in combined_dl_src.keys():
combined_dl_src[results['name'][isource]] = [name, exclArea]
else:
if combined_dl_src[results['name'][isource]][1] > exclArea:
combined_dl_src[results['name'][isource]] = [name, exclArea]
# make sure you have removed all
# unwanted sources from likelihood combination
for n in no_lim:
combined_dl_src.pop(n, None)
minArea, maxArea = np.min(exclAreas.values()),np.max(exclAreas.values())
m = 1. / (maxArea - minArea)
b = -minArea * m
In [19]:
plt.close()
fig = plt.figure(0)
ax = fig.add_subplot(111)
handles = []
for isource, name in enumerate(results['SOURCE_FULL']):
LL, BB = np.meshgrid(10.**results['lcoh'][isource][:,0],
10.**results['igmf'][isource][0,:],
indexing = 'ij')
dloglike = results['dloglike'][isource]
if results['name'][isource] in no_lim: continue
plt.contourf(LL,BB,2*dloglike, [-1e5,-6.18],
colors = [cp( m * exclAreas[name] + b)],
linestyles = ['-'],
zorder = -(m * exclAreas[name] + b),
alpha = 0.5
)
c = plt.contour(LL,BB,2*dloglike, [-1e5,-6.18],
colors = [cp( m * exclAreas[name] + b - 0.05)],
linestyles = ['--'],
label = name
)
handles.append(mpatches.Patch(
color= cp( m * exclAreas[name] + b),
alpha = 0.5,
label = ' '.join(name.split('_'))))
title = '$t_\mathrm{{max}} = 10^{0:.0f}$ yrs, $\\theta_\mathrm{{jet}} = {1:.1f}^\circ$'.format(
np.log10(config['tmax']), config['th_jet'])
plt.grid(True)
plt.legend(loc = 'lower left', handles = handles, fontsize = 'x-small',
ncol = 2, title = title)
ax.get_legend().get_title().set_fontsize('small')
ax.set_xscale('log')
ax.set_yscale('log')
ax.set_ylabel('$B\,$(G)')
ax.set_xlabel('$\lambda\,$(Mpc)')
fname = 'dL_tmax{0[tmax]:.0e}_thj{0[th_jet]:.0f}_{0[kind]:s}'.format(config)
plt.savefig('./plots/' + fname + '.png', format = 'png', dpi = 200)
plt.savefig('./plots/' + fname + '.pdf', format = 'pdf')
In [20]:
plt.close()
fig = plt.figure(1)
ax = fig.add_subplot(111)
handles = []
for isource, name in enumerate(results['SOURCE_FULL']):
LL, BB = np.meshgrid(10.**results['lcoh'][isource][:,0],
10.**results['igmf'][isource][0,:],
indexing = 'ij')
dloglike = np.sum(results['loglike_comp'][isource][...,1:], axis = -1)
dloglike = -(dloglike - np.min(dloglike))
if results['name'][isource] in no_lim: continue
plt.contourf(LL,BB,2*dloglike, [-1e10,-6.18],
colors = [cp( m * exclAreas[name] + b)],
linestyles = ['-'],
zorder = -(m * exclAreas[name] + b),
alpha = 0.5
)
c = plt.contour(LL,BB,2*dloglike, [-1e10,-6.18],
colors = [cp( m * exclAreas[name] + b - 0.05)],
linestyles = ['--'],
label = name
)
handles.append(mpatches.Patch(
color= cp( m * exclAreas[name] + b),
alpha = 0.5,
label = ' '.join(name.split('_'))))
title = '$t_\mathrm{{max}} = 10^{0:.0f}$ yrs, $\\theta_\mathrm{{jet}} = {1:.1f}^\circ$'.format(
np.log10(config['tmax']), config['th_jet'])
plt.grid(True)
plt.legend(loc = 'lower left', handles = handles, fontsize = 'x-small',
ncol = 2, title = title)
ax.get_legend().get_title().set_fontsize('small')
ax.set_xscale('log')
ax.set_yscale('log')
ax.set_ylabel('$B\,$(G)')
ax.set_xlabel('$\lambda\,$(Mpc)')
fname = 'dLFermi_tmax{0[tmax]:.0e}_thj{0[th_jet]:.0f}_{0[kind]:s}'.format(config)
plt.savefig('./plots/' + fname + '.png', format = 'png', dpi = 200)
plt.savefig('./plots/' + fname + '.pdf', format = 'pdf')
In [21]:
combined_dL_spec = [x[0] for x in combined_dl_src.values()]
dloglike = []
for isource, name in enumerate(results['SOURCE_FULL']):
if not name in combined_dL_spec: continue
dloglike.append(results['dloglike'][isource])
dloglike = np.array(dloglike).sum(axis = 0)
Plot the combined likelihood:
In [22]:
ax = plt.subplot(111)
plt.contourf(LL,BB,2*dloglike, [-1e5,-6.18],
colors = ['0.5'],
linestyles = ['-'],
alpha = 0.5
)
title = '$t_\mathrm{{max}} = 10^{0:.0f}$ yrs, $\\theta_\mathrm{{jet}} = {1:.1f}^\circ$'.format(
np.log10(config['tmax']), config['th_jet'])
title += '\nSources:'
for c in combined_dL_spec:
c = c.split('_')
title += '\n{0[0]:s} {0[1]:s}'.format(c)
handles = [mpatches.Patch(
color= '0.5',
alpha = 0.5,
label = 'Combined 95% Exclusion')]
plt.legend(loc = 4, title = title, handles = handles)
plt.grid(True)
ax.set_xscale('log')
ax.set_yscale('log')
ax.set_ylabel('$B\,$(G)')
ax.set_xlabel('$\lambda\,$(Mpc)')
fname = 'Comb_dL_tmax{0[tmax]:.0e}_thj{0[th_jet]:.0f}_{0[kind]:s}'.format(config)
plt.savefig('./plots/' + fname + '.png', format = 'png', dpi = 200)
plt.savefig('./plots/' + fname + '.pdf', format = 'pdf')
In [23]:
dloglike = results['dloglike'][isource]
#plt.pcolormesh(10**ll,10**bb,2*dloglike,
# cmap=cmap, #edgecolors = 'skyblue', lw = 0.01,
# vmax=0.0,vmin=int(np.min(dloglike))
# )
#cb = plt.colorbar(label = '$2\Delta \log \mathcal{L}$')
#plt.gca().set_xlim(np.min(10.**ll),np.max(10.**ll))
#plt.gca().set_ylim(np.min(10.**bb),np.max(10.**bb))
cp = plt.cm.gnuplot
cp.set_over('1.')
BB, LL = np.meshgrid(10.**results['igmf'][isource][0,:], 10.**results['lcoh'][isource][:,0])
plot_map_countours(LL,BB,2*dloglike, [-11.83,-6.18],
colbarlabel = '$2\Delta \log \mathcal{L}$',
colmap = cp,
#vmin = -50, vmax = 0.,
vmin = np.min(2*dloglike), vmax = -6.
#set_under = '1.'
)
plt.gca().set_xscale('log')
plt.gca().set_yscale('log')
plt.gca().set_ylabel('$B\,$(G)')
plt.gca().set_xlabel('$\lambda\,$(Mpc)')
tics_out(plt.gca())
Out[23]:
In [24]:
loglike = np.rollaxis(results['loglike_comp'][isource],2)
label = ['IACT SED', '$\mathit{Fermi}$ SED', 'Cascade']
plt.figure(figsize = (6,4*3))
for i,logl in enumerate(loglike):
ax = plt.subplot(loglike.shape[0],1,i+1)
dl = -2.*(logl - np.min(logl))
#dl = -logl
plt.pcolormesh(10**ll,10**bb,
dl,
cmap=cp,
# edgecolors = 'skyblue', lw = 0.01,
vmax=int(np.max(dl)),
vmin=int(np.min(dl)))
cb = plt.colorbar(label = '$2\Delta \log \mathcal{L}$')
ax.set_ylabel('$B\,$(G)')
ax.set_xlabel('$\lambda\,$(Mpc)')
#ax.grid(True, color = 'skyblue')
ax.set_xscale('log')
ax.set_yscale('log')
tics_out(ax)
ax.set_xlim(np.min(10.**ll),np.max(10.**ll))
ax.set_ylim(np.min(10.**bb),np.max(10.**bb))
ax.annotate(label[i], xy = (0.05,0.05), xycoords = 'axes fraction', **effect)
plt.subplots_adjust(hspace = 0.3)
In [25]:
dloglike_fermi = np.sum(results['loglike_comp'][isource][...,1:], axis = -1)
#dloglike -= np.max(dloglike)
dloglike_fermi = -(dloglike_fermi - np.min(dloglike_fermi))
plot_map_countours(LL,BB,2*dloglike_fermi, [-11.83,-6.18],
colbarlabel = '$2\Delta \log \mathcal{L}$',
colmap = cp,
#vmin = -50, vmax = 0.,
vmin = np.min(2*dloglike_fermi), vmax = -6.
#set_under = '1.'
)
plt.gca().set_xscale('log')
plt.gca().set_yscale('log')
plt.gca().set_ylabel('$B\,$(G)')
plt.gca().set_xlabel('$\lambda\,$(Mpc)')
#plt.gca().grid(True, color = 'skyblue')
#plt.gca().set_xlim(np.min(10.**ll),np.max(10.**ll))
#plt.gca().set_ylim(np.min(10.**bb),np.max(10.**bb))
tics_out(plt.gca())
Out[25]:
In [ ]: