<jupyter_start><jupyter_text>Notebook to test the Fit of SEDs to an IGMF model Inspired from https://github.com/woodmd/haloanalysis/blob/master/haloanalysis/scripts/fit_igmf.pyL65-L66<jupyter_code>%matplotlib inline
%load_ext autoreload
%autoreload 2<jupyter_output><empty_output><jupyter_text>Imports<jupyter_code>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<jupyter_output><empty_output><jupyter_text>load the catalog Catalog file:<jupyter_code>cat = '../data/table_std_psf0123_joint2a_stdmodel_lnl_v14.fits'
tab_pars = Table.read(cat,hdu='SCAN_PARS')
tab_ebounds = Table.read(cat,hdu='EBOUNDS')
tables = [Table.read(cat,'CATALOG'),
            Table.read(cat,'LIKELIHOOD'),
            Table.read(cat,'SED')]
for i, t in enumerate(tables):
    if 'NAME' in t.columns:
        t['name'] = t['NAME']
print tables[0]['name'].data
print tables[1]['name'].data
print tables[2]['name'].data<jupyter_output>['FHES J2359.8-3736' 'FHES J0001.2-0746' 'FHES J0001.5+2113' ...,
 'FHES J0538.0-6913' 'FHES J0535.6-6734' 'FHES J0524.7-6937']
['3FGL J0000.2-3738' '3FGL J0001.2-0748' '3FGL J0001.4+2120' ..., 'LMC P2'
 'LMC P3' 'LMC P4']
['3FGL J0000.2-3738' '3FGL J0001.2-0748' '3FGL J0001.4+2120' ..., 'LMC P2'
 'LMC P3' 'LMC P4']<jupyter_text>Remove identifier so in names so that join operation works<jupyter_code>for i, t in enumerate(tables):
    for j,n in enumerate(t['name']):
        tables[i]['name'][j] = n[5:-2] ### WORK AROUND
print tables[0]['name'].data<jupyter_output>['J2359.8-37' 'J0001.2-07' 'J0001.5+21' ..., 'J0538.0-69' 'J0535.6-67'
 'J0524.7-69']<jupyter_text>join the tables -- DOES ONLY WORK WITH WORKAROUND ABOVE<jupyter_code>tab_casc = join(tables[0],tables[1]) # table is empty since 'name' names are different
tab_casc = join(tab_casc,tables[2])
print tab_casc['name'].data.shape<jupyter_output>(638,)<jupyter_text>Load the TeV SEDs<jupyter_code>tab_sed_tev = Table.read('../data/CompiledTeVSources.fits')
print tab_sed_tev<jupyter_output>3FGL_NAME             SOURCE_FULL            ...   NORM_ERRN [25]  
------------ -------------------------------- ... ------------------
J0232.8+2016       1ES0229+200_HESS_2005-2006 ...  1.3277e-08 .. nan
J0232.8+2016    1ES0229+200_VERITAS_2009-2012 ...     2.8e-08 .. nan
J0349.2-1158      1ES0347-121_HESS_2006-09-12 ...  3.5758e-08 .. nan
J0416.8+0104       1ES0414+009_HESS_2005-2009 ...   1.017e-07 .. nan
J0416.8+0104    1ES0414+009_VERITAS_2008-2011 ...    1.09e-08 .. nan
J0809.8+5218    1ES0806+524_VERITAS_2006-2008 ... 9.29454e-08 .. nan
J1015.0+4925           1ES1011+496_MAGIC_2007 ...    1.69e-06 .. nan
J1103.5-2329       1ES1101-232_HESS_2004-2005 ...  1.0851e-07 .. nan
J1217.8+3007           1ES1215+303_MAGIC_2011 ...    1.53e-06 .. nan
J1217.8+3007    1ES1215+303_VERITAS_2008-2012 ...  5.7133e-08 .. nan
         ...                              ... ...                ...
J2158.8-3013      PKS2155-304_HESS_2006-08-03 ...   2.947e-07 .. nan
J2158.8-3013      PKS2155-304_HESS[...]<jupyter_text>Load the IGMF model cube<jupyter_code>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',
)
casc_model.set_eblmodel(eblmodel = 'dominguez')<jupyter_output><empty_output><jupyter_text>Loop over the sources and perform the likelihood fit Spectra you want to include in the analysis<jupyter_code>src_names = ['1ES0229+200_HESS_2005-2006']
nstep = 5
casc_scale = 1.
casc_r68_scale = 1.
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]<jupyter_output>3FGL_NAME          SOURCE_FULL         ...   NORM_ERRP [25]    NORM_ERRN [25] 
------------ -------------------------- ... ----------------- -----------------
J0232.8+2016 1ES0229+200_HESS_2005-2006 ... 1.3277e-08 .. nan 1.3277e-08 .. nan
   name         codename     ...       dloglike_scan [20,9]     
                             ...                                
---------- ----------------- ... -------------------------------
J0232.8+20 3fgl_j0232.8+2016 ... -4.2029670149 .. -2.70594689762
   name         codename     ...       dloglike_scan [20,9]     
                             ...                                
---------- ----------------- ... -------------------------------
J0232.8+20 3fgl_j0232.8+2016 ... -4.2029670149 .. -2.70594689762
((0, 0), -65.160357820170816, [0.14000000000000001, -4.0, -20.0], [1.7781826364287957e-13, -1.6174113205437628, 4929568.6279378785])
((0, 1), -66.836542031590739, [0.14000000000000001, -4.0, -18.0], [1.7762412008607477e-13, -1.6170279304050[...]
{'fraction_non_alphanumeric': 0.17637822749476623, 'fraction_numerical': 0.15247732030704816, 'mean_word_length': 2.5853658536585367, 'pattern_counts': {'":': 0, '<': 28, '<?xml version=': 0, '>': 28, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 2, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '600773', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 2748, 'n_tokens_neox': 2241, 'n_words': 375, 'path': 'notebooks/fit_igmf_test.ipynb', 'repo_name': 'fermiPy/haloanalysis', 'size': '5732'}