In [1]:
%matplotlib inline
In [2]:
%load_ext autoreload
In [3]:
%autoreload 2
In [4]:
import os
import sys
from collections import OrderedDict
import yaml
import numpy as np
from astropy.io import fits
from astropy.table import Table, Column, join, hstack, vstack
import matplotlib.pyplot as plt
import matplotlib
import fermipy.utils as utils
from fermipy.spectrum import *
from fermipy.castro import CastroData
from haloanalysis.utils import create_mask, load_source_rows
from haloanalysis.sed import HaloSED
from haloanalysis.model import CascModel, CascLike
from haloanalysis.model import scan_igmf_likelihood
Catalog file:
In [5]:
cat = '../data/table_std_psf0123_joint2a_stdmodel_lnl_v14.fits'
In [6]:
tab_pars = Table.read(cat,hdu='SCAN_PARS')
tab_ebounds = Table.read(cat,hdu='EBOUNDS')
In [7]:
tables = [Table.read(cat,'CATALOG'),
Table.read(cat,'LIKELIHOOD'),
Table.read(cat,'SED')]
In [8]:
for i, t in enumerate(tables):
if 'NAME' in t.columns:
t['name'] = t['NAME']
In [9]:
print tables[0]['name'].data
print tables[1]['name'].data
print tables[2]['name'].data
Remove identifier so in names so that join operation works
In [10]:
for i, t in enumerate(tables):
for j,n in enumerate(t['name']):
tables[i]['name'][j] = n[5:-2] ### WORK AROUND
In [11]:
print tables[0]['name'].data
In [12]:
tab_casc = join(tables[0],tables[1]) # table is empty since 'name' names are different
tab_casc = join(tab_casc,tables[2])
In [13]:
print tab_casc['name'].data.shape
Load the TeV SEDs
In [14]:
tab_sed_tev = Table.read('../data/CompiledTeVSources.fits')
In [15]:
print tab_sed_tev
In [17]:
casc_model = CascModel.create_from_fits(
'/nfs/farm/g/glast/u/mmeyer/projects/FermiHalo/Output/EBLm6/th_jet6.00/gam-2.00/results_merged_z_th6d_t1e7.fits',
)
In [20]:
casc_model.set_eblmodel(eblmodel = 'dominguez')
Spectra you want to include in the analysis
In [21]:
src_names = ['1ES0229+200_HESS_2005-2006']
In [22]:
nstep = 5
casc_scale = 1.
casc_r68_scale = 1.
In [23]:
tab_igmf = []
for name in src_names:
rows_sed_tev = load_source_rows(tab_sed_tev, [name], key='SOURCE_FULL')
#cat_names = [ 'FHES %s'%row['3FGL_NAME'] for row in rows_sed_tev ]
cat_names = [ '%s'%row['3FGL_NAME'][:-2] for row in rows_sed_tev ] # WORK AROUND
cat_names = np.unique(np.array(cat_names))
rows_sed_gev = load_source_rows(tab_casc, cat_names, key='name')
rows_casc = load_source_rows(tab_casc, cat_names, key='name')
print rows_sed_tev
print rows_sed_gev
print rows_casc
tab = scan_igmf_likelihood(casc_model, rows_sed_tev, rows_sed_gev,
rows_casc, tab_pars, tab_ebounds, nstep,
casc_scale, casc_r68_scale)
#tab_igmf += [tab]
((0, 0), -65.160357820170816, [0.14000000000000001, -4.0, -20.0], [1.7781826364287957e-13, -1.6174113205437628, 4929568.6279378785])
In [ ]: