In [3]:

    

import pdb
import numpy as np
import pandas as pd
import os
import math
import matplotlib.pyplot as plt
import pylab as plt
from astropy.wcs import WCS
from astropy.io import fits
import readcol
from utils import fast_sed_fitter
from utils import fast_double_sed_fitter
from utils import fast_Lir
from utils import fast_double_Lir
from utils import L_fun
from utils import L_fit
from utils import T_fun
from utils import T_fit
from utils import shift_twod
from utils import dist_idl
from utils import loggen
from utils import gauss_kern
from lmfit import Parameters, minimize, fit_report
from simstack import stack_libraries_in_redshift_slices
from skymaps import Skymaps
from skymaps import Field_catalogs

%matplotlib inline


    

In [4]:

    

popcolor=['blue','red','green','orange','black','grey','chocolate','darkviolet','pink',
          'magenta','dodgerblue','lavender','blue','red','green','orange','black','grey',
          'chocolate','darkviolet','pink','magenta','dodgerblue','lavender']


    

In [5]:

    

conv_sfr = 1.728e-10 / 10**(.23)


    

In [6]:

    

#Define Redshift Bins
z_lo = np.array([0.0,0.5,1.0,1.5,2.0,2.5,3.0,3.5])
z_hi = np.array([0.5,1.0,1.5,2.0,2.5,3.0,3.5,4.0])
z_nodes = np.array([0.5,1.0,1.5,2.0,2.5,3.0,3.5,4.0])
z_mid = (z_nodes[:-1] + z_nodes[1:])/2
nz = len(z_nodes) - 1


    

In [7]:

    

#Define Stellar Mass Bins
m_lo = np.array([8.5, 9.5,10.0,10.5,11.0])
m_hi = np.array([9.5,10.0,10.5,11.0,13.0])
m_nodes = np.array([8.5, 9.5,10.0,10.5,11.0,13.0])
m_mid = (m_nodes[:-1] + m_nodes[1:])/2
nm = len(m_nodes) - 1


    

In [8]:

    

#Decide which Maps to include
wv0 = np.array([1,1,1,1,1,1,1,1,1]) # This includes all maps
wv0 = np.array([0,0,1,1,1,0,1,1,0]) # 
wv0 = np.array([0,0,1,1,1,0,1,1,1]) # 
indstack = np.where(wv0 == 1)


    

In [9]:

    

## Map Directories 

dir_spitzer_maps = '/home/agn/viero/Jason/cosmos/'
dir_spire_maps = '/home/agn/viero/Jason/cosmos/'
dir_pacs_maps = '/home/agn/viero/Jason/cosmos/'
dir_scuba_maps = '/home/agn/viero/Jason/cosmos/'
dir_aztec_maps = '/home/agn/viero/Jason/cosmos/'
dir_aztec_maps = '/home/agn/viero/Jason/cosmos/'

## Dictionary Names
library_keys =['mips24'
               ,'pacs_green'
               ,'pacs_red'
               ,'spire_PSW'
               ,'spire_PMW'
               ,'scuba_450'
               ,'spire_PLW'
               ,'scuba_850'
               ,'aztec'
              ]

wavelength=[24.,100,160,250,350,450,500,850,1100]
nwv = np.sum(wv0) 
fwhm =[6.32,7.4, 11.3,18.1, 25.2,7., 36.6,15.,18.]
efwhm=[6.32,6.7, 11.2,17.6, 23.9,7.8, 35.2,14.5,18.] # want to the measured effective FWHM later
color_correction=[1.25,23.58,23.82,1.018,0.9914,1e-3,0.95615,1e-3,1.0]
beam_area = [1.55e-09,1.,1.,1.,1.,1.62e-09,1.,5.6e-09,1.] #sr
beam_area = [1.55e-09,1.,1.,1.,1.,1.,1.,1.,1.] #sr

pixsize_suffix = '6.0_arcsec_pixels'

maps = [dir_spitzer_maps+'mips_24_GO3_sci_10.cutout.fits' 
        ,dir_pacs_maps+'pep_COSMOS_green_Map.DR1.sci.cutout.fits'
        ,dir_pacs_maps+'pep_COSMOS_red_Map.DR1.sci.cutout.fits'
        ,dir_spire_maps+'cosmos-uvista-hipe12_itermap_10_iterations_'+pixsize_suffix+'_PSW.signal.cutout.fits'
        ,dir_spire_maps+'cosmos-uvista-hipe12_itermap_10_iterations_'+pixsize_suffix+'_PMW.signal.cutout.fits'
        ,dir_scuba_maps+'map450_new_header.cutout.fits'
        ,dir_spire_maps+'cosmos-uvista-hipe12_itermap_10_iterations_'+pixsize_suffix+'_PLW.signal.cutout.fits'
        ,dir_scuba_maps+'map850_new_header.cutout.fits'
        ,dir_aztec_maps+'cosmos_jcmt_kscott20100925_new_header_6.0_arcsec_pixels_map.fits'
       ]
noises = [dir_spitzer_maps+'mips_24_GO3_unc_10.cutout.fits' 
        ,dir_pacs_maps+'pep_COSMOS_green_Map.DR1.err.cutout.fits'
        ,dir_pacs_maps+'pep_COSMOS_red_Map.DR1.err.cutout.fits'
        ,dir_spire_maps+'cosmos-uvista-hipe12_itermap_10_iterations_'+pixsize_suffix+'_PSW.noise.cutout.fits'
        ,dir_spire_maps+'cosmos-uvista-hipe12_itermap_10_iterations_'+pixsize_suffix+'_PMW.noise.cutout.fits'
        ,dir_scuba_maps+'map450_new_header_rms.cutout.fits'
        ,dir_spire_maps+'cosmos-uvista-hipe12_itermap_10_iterations_'+pixsize_suffix+'_PLW.noise.cutout.fits'
        ,dir_scuba_maps+'map850_new_header_rms.cutout.fits'
        ,dir_aztec_maps+'cosmos_jcmt_kscott20100925_new_header_6.0_arcsec_pixels_noise.fits'
       ]


    

In [ ]:

    




    

In [10]:

    

sky_library={}

for t in indstack[0]:
    sky_library[library_keys[t]] = Skymaps(maps[t],noises[t],efwhm[t],color_correction=color_correction[t])
    sky_library[library_keys[t]].add_wavelength(wavelength[t])
    sky_library[library_keys[t]].add_fwhm(efwhm[t]) 
    if beam_area[t] != 1: sky_library[library_keys[t]].beam_area_correction(beam_area[t])


    

    




---------------------------------------------------------------------------
IOError                                   Traceback (most recent call last)
<ipython-input-10-c882a52f6e37> in <module>()
      2 
      3 for t in indstack[0]:
----> 4     sky_library[library_keys[t]] = Skymaps(maps[t],noises[t],efwhm[t],color_correction=color_correction[t])
      5     sky_library[library_keys[t]].add_wavelength(wavelength[t])
      6     sky_library[library_keys[t]].add_fwhm(efwhm[t])

/home/marcoviero/Code/Python/Modules/Utils/skymaps.py in __init__(self, file_map, file_noise, psf, color_correction)
     31                 else:
     32                         #This assumes that if Signal and Noise are different maps, they are contained in first extension
---> 33                         cmap, hd = fits.getdata(file_map, 0, header = True)
     34                         cnoise, nhd = fits.getdata(file_noise, 0, header = True)
     35 

/home/marcoviero/miniconda2/lib/python2.7/site-packages/astropy/io/fits/convenience.pyc in getdata(filename, *args, **kwargs)
    199     view = kwargs.pop('view', None)
    200 
--> 201     hdulist, extidx = _getext(filename, mode, *args, **kwargs)
    202     try:
    203         hdu = hdulist[extidx]

/home/marcoviero/miniconda2/lib/python2.7/site-packages/astropy/io/fits/convenience.pyc in _getext(filename, mode, *args, **kwargs)
    887         raise TypeError('extver alone cannot specify an extension.')
    888 
--> 889     hdulist = fitsopen(filename, mode=mode, **kwargs)
    890 
    891     return hdulist, ext

/home/marcoviero/miniconda2/lib/python2.7/site-packages/astropy/io/fits/hdu/hdulist.pyc in fitsopen(name, mode, memmap, save_backup, cache, lazy_load_hdus, **kwargs)
    164 
    165     return HDUList.fromfile(name, mode, memmap, save_backup, cache,
--> 166                             lazy_load_hdus, **kwargs)
    167 
    168 

/home/marcoviero/miniconda2/lib/python2.7/site-packages/astropy/io/fits/hdu/hdulist.pyc in fromfile(cls, fileobj, mode, memmap, save_backup, cache, lazy_load_hdus, **kwargs)
    402         return cls._readfrom(fileobj=fileobj, mode=mode, memmap=memmap,
    403                              save_backup=save_backup, cache=cache,
--> 404                              lazy_load_hdus=lazy_load_hdus, **kwargs)
    405 
    406     @classmethod

/home/marcoviero/miniconda2/lib/python2.7/site-packages/astropy/io/fits/hdu/hdulist.pyc in _readfrom(cls, fileobj, data, mode, memmap, save_backup, cache, lazy_load_hdus, **kwargs)
   1013             if not isinstance(fileobj, _File):
   1014                 # instantiate a FITS file object (ffo)
-> 1015                 fileobj = _File(fileobj, mode=mode, memmap=memmap, cache=cache)
   1016             # The pyfits mode is determined by the _File initializer if the
   1017             # supplied mode was None

/home/marcoviero/miniconda2/lib/python2.7/site-packages/astropy/utils/decorators.pyc in wrapper(*args, **kwargs)
    491                         # one with the name of the new argument to the function
    492                         kwargs[new_name[i]] = value
--> 493             return function(*args, **kwargs)
    494 
    495         return wrapper

/home/marcoviero/miniconda2/lib/python2.7/site-packages/astropy/io/fits/file.pyc in __init__(self, fileobj, mode, memmap, overwrite, cache)
    142             self._open_fileobj(fileobj, mode, overwrite)
    143         elif isinstance(fileobj, string_types):
--> 144             self._open_filename(fileobj, mode, overwrite)
    145         else:
    146             self._open_filelike(fileobj, mode, overwrite)

/home/marcoviero/miniconda2/lib/python2.7/site-packages/astropy/io/fits/file.pyc in _open_filename(self, filename, mode, overwrite)
    497             self._file = bz2.BZ2File(self.name, bzip2_mode)
    498         else:
--> 499             self._file = fileobj_open(self.name, IO_FITS_MODES[mode])
    500 
    501         # Make certain we're back at the beginning of the file

/home/marcoviero/miniconda2/lib/python2.7/site-packages/astropy/io/fits/util.pyc in fileobj_open(filename, mode)
    391         """
    392 
--> 393         return open(filename, mode)
    394 else:
    395     def fileobj_open(filename, mode):

IOError: [Errno 2] No such file or directory: '/home/agn/viero/Jason/cosmos/pep_COSMOS_red_Map.DR1.sci.cutout.fits'

In [ ]:

    




    

In [ ]:

    




    

In [10]:

    

path_catalog = '/data/maps_cats_models/catalogs/UVISTA/'
#Catalog must have the following (case sensitive) Columns:
# - RA
# - DEC
# - z_peak
# - LMASS
# - rf_U_V
# - rf_V_J
#ALPHA_J2000,DELTA_J2000,ZPDF,ZPDF_L68,ZPDF_H68,CLASS,MASS_MED,MASS_MED_MIN68,MASS_MED_MAX68
file_catalog = 'COSMOS2015_Laigle+_v1.1.csv'


    

In [11]:

    

tbl = pd.read_table(path_catalog+file_catalog,sep=',')
tbl['ID']=range(len(tbl['ALPHA_J2000']))
tbl['ra']=tbl['ALPHA_J2000']
tbl['dec']=tbl['DELTA_J2000']
tbl['z_peak']=tbl['ZPDF']
tbl['sfg']=tbl['CLASS']
tbl['LMASS']=tbl['MASS_MED']


    

In [12]:

    

tbl.keys()


    

    
Out[12]:





Index([u'ALPHA_J2000', u'DELTA_J2000', u'ZPDF', u'ZPDF_L68', u'ZPDF_H68',
       u'CLASS', u'MASS_MED', u'MASS_MED_MIN68', u'MASS_MED_MAX68', u'ID',
       u'ra', u'dec', u'z_peak', u'sfg', u'LMASS'],
      dtype='object')

In [13]:

    

uVista = Field_catalogs(tbl)


    

In [ ]:

    

#Stack in redshift bins, with Layers divided by Stellar Mass and Star-Forming/Quiescent Galaxies
pop = ['sf','qt']
npop=2
all_stacked_fluxes = np.zeros([nwv,nz,nm,npop])
all_luminosity_temp = np.zeros([nz,nm,npop,2])

for iz in range(nz):
    zn = z_nodes[iz:iz+2]
    uVista.get_sf_qt_mass_redshift_bins(zn,m_nodes)
    radec_m_z_p = uVista.subset_positions(uVista.id_z_ms)
    stacked_fluxes_  =  None
    n_sources_max = None
    
    stacked_fluxes = stack_libraries_in_redshift_slices(
        sky_library,
        radec_m_z_p)
    
    #pdb.set_trace()
    args = radec_m_z_p.keys()
    for iwv in range(nwv):
        stacked_fluxes_wv = stacked_fluxes[str(wavelength[indstack[0][iwv]])]
        for j in range(nm):
            for k in range(npop):
                arg = 'z_'+str(zn[0])+'-'+str(zn[1])+'__m_'+str(m_nodes[j])+'-'+str(m_nodes[j+1])+'_'+pop[k]
                all_stacked_fluxes[iwv,iz,j,k] = stacked_fluxes_wv[arg.replace('.','p').replace('-','_')].value
    
    #pdb.set_trace()
    #PLOT
    plt.figure()
    plt.ylim([-1e-1,3e0])
    plt.ylim([1e-2,9e1])
    plt.xlim([20,510])
    plt.xlim([20,1800])
    plt.xlim([20,2500])
    plt.yscale('log')
    plt.ylabel('Flux Density [mJy]',fontsize=14)
    plt.xlabel('Wavelength [microns]',fontsize=14)
    ln=['','--']
    for k in range(npop):
        for j in range(nm)[::-1]: 
            arg = 'z_'+str(round(zn[0],3))+'-'+str(round(zn[1],3))+'__m_'+str(round(m_nodes[j],3))+'-'+str(round(m_nodes[j+1],3))+'_'+pop[k]
            ng = len(radec_m_z_p[arg.replace('.','p').replace('-','_')][0])
            plt.plot(np.array(wavelength)[indstack],1e3*all_stacked_fluxes[:,iz,j,k],ln[k],color=popcolor[npop*j+k],label=[str(ng)+' '+pop[k],str(m_mid[j])])
            plt.legend()
    plt.show()
    #pdb.set_trace()


    

In [ ]:

    

def nghack(a):
    b = a
    for i in range(len(a)):
        if a[i] < 0:
            b[i] = 1e-5
    return b


    

In [ ]:

    

for iwv in range(nwv):
    plt.figure()
    for j in range(nm):
        plt.yscale('log')
        plt.ylabel('Flux Density [mJy]',fontsize=14)
        plt.xlabel('Redshift',fontsize=14)
        plt.title('Star-Forming Galaxies, '+str(wavelength[indstack[0][iwv]]))
        plt.plot(z_mid,1e3*nghack(all_stacked_fluxes[iwv,:,j,0]),'o')
        plt.plot(z_mid,1e3*nghack(all_stacked_fluxes[iwv,:,j,0]))
    plt.show()


    

In [ ]:

    

for iwv in range(nwv):
    plt.figure()
    for j in range(nm):
        plt.yscale('log')
        plt.ylabel('Flux Density [mJy]',fontsize=14)
        plt.xlabel('Redshift',fontsize=14)
        plt.title('Star-Forming Galaxies, '+str(wavelength[indstack[0][iwv]]))
        plt.plot(z_mid,1e3*nghack(all_stacked_fluxes[iwv,:,j,0]),'o')
        plt.plot(z_mid,1e3*nghack(all_stacked_fluxes[iwv,:,j,0]))
    plt.show()


    

In [ ]:

    

for iwv in range(nwv):
    plt.figure()
    for j in range(nm):
        plt.yscale('log')
        plt.ylim([1e-2,9e1])
        plt.ylabel('Flux Density [mJy]',fontsize=14)
        plt.xlabel('Redshift',fontsize=14)
        plt.title('Quiescent Galaxies, '+str(wavelength[indstack[0][iwv]]))
        plt.plot(z_mid,1e3*(all_stacked_fluxes[iwv,:,j,1]),'o')
        plt.plot(z_mid,1e3*(all_stacked_fluxes[iwv,:,j,1]))
    plt.show()


    

In [ ]:

    

wvs = np.array(wavelength)[indstack]
luminosities = np.zeros([nz,nm])
temperatures = np.zeros([nz,nm])
luminosity_dictionary = {}
lumfit_dictionary = {}

for iz in range(nz):
    zn = z_nodes[iz:iz+2] 
    plt.ylim([9,14])
    for j in range(nm)[::-1]:
        flux_wv = (all_stacked_fluxes[:,iz,j,0] )   
        m = fast_sed_fitter(wvs,flux_wv,covar=np.sqrt(abs(flux_wv)))
        #sed = fast_sed(m, wvs)
        L_mod = fast_Lir(m, np.mean(zn))
        luminosities[iz,j] = L_mod.value
        temperatures[iz,j] = m['T_observed'].value * np.mean(zn)
        #print temperatures[iz,k] 
        
for j in range(nm)[::-1]:
    luminosity_dictionary['L_sed'+str(j)] = luminosities[:,j]
    luminosity_dictionary['T_sed'+str(j)] = temperatures[:,j]
        #print LIR.value
        #plt.plot(np.mean(zn),np.log10(LIR.value),'o')
    plt.plot(z_mid,np.log10(luminosity_dictionary['L_sed'+str(j)]))
        #plt.yscale('log')
        
    plt.ylabel('Luminosity [L_star]',fontsize=14)
    plt.xlabel('Redshift',fontsize=14)
#plt.show()


    

In [ ]: