# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python (herschelhelp_internal)
#     language: python
#     name: helpint
# ---

# This notebook prepare the catalogues that will be analysed by CIGALE for SED fitting and physical parameter estimation.

# +
import numpy as np
import os
os.environ['LOG_LEVEL'] = 'INFO'
import glob
from astropy.table import Table

from herschelhelp.filters import correct_galactic_extinction
from herschelhelp.external import convert_table_for_cigale
# -

SUFFIX = '20180221'
file_list = glob.glob('../../dmu32/dmu32_SGP/data/tiles/sgp_{}_cigale_*.fits'.format(SUFFIX))

# +
#len(master_catalogue)
# -

# # Best sources
#
# Define a good far-IR measurement as:
# - an existing flux in the band;
# - the flag from XID+ must not be set;
# - the signal to noise ratio must be over 2.

# +
# good = {}
# for band in ['pacs_green', 'pacs_red', 'spire_250', 'spire_350', 'spire_500']:
#     good[band] = (~np.isnan(master_catalogue['f_{}'.format(band)]) & 
#                                    ~master_catalogue['flag_{}'.format(band)])
#     good[band][good[band]] &= (master_catalogue[good[band]]['f_{}'.format(band)] /
#                               master_catalogue[good[band]]['ferr_{}'.format(band)] >= 2)
# -

# We will keep only sources with at leat 2 good far-IR measurements (we may actually use less sources are not all may have a redshift).

# +
# combined_good = np.sum(list(good.values()), axis=0) >= 2

# +
# print("Number of good sources: {}".format(np.sum(combined_good)))
# -

# # Main catalogue for CIGALE

# +
# best_catalogue = master_catalogue[combined_good].copy()

# +
# # Correction for galactic extinction
# best_catalogue = correct_galactic_extinction(best_catalogue, inplace=True)

# +
# # Convertion to CIGALE format
# best_catalogue = convert_table_for_cigale(best_catalogue, inplace=True, remove_zerofluxes=True)
# -

# ## Band selection
#
# We want to use only one filter for similar bands. We define an order of preference and set to NaN the flux in the lower prefered bands when a prefered band is available. Some band may have a 0 flux, we set there values to NaN.

# +
u_bands = ["omegacam_u"]
g_bands = ["omegacam_g", "decam_g",            "gpc1_g"]
r_bands = ["omegacam_r", "decam_r",            "gpc1_r"]
i_bands = ["omegacam_i", "decam_i",            "gpc1_i"]
z_bands = [              "decam_z", "vista_z", "gpc1_z"]
y_bands = [              "decam_y", "vista_y", "gpc1_y"]

def remove_unneeded_fluxes(list_of_bands):
    for band_idx, band in enumerate(list_of_bands[:-1]):
        mask = ~np.isnan(best_catalogue[band])
        for lower_band in list_of_bands[band_idx+1:]:
            best_catalogue[lower_band][mask] = np.nan
            best_catalogue["{}_err".format(lower_band)][mask] = np.nan


# +
# remove_unneeded_fluxes(g_bands)
# remove_unneeded_fluxes(u_bands)
# remove_unneeded_fluxes(r_bands)
# remove_unneeded_fluxes(i_bands)
# remove_unneeded_fluxes(z_bands)
# remove_unneeded_fluxes(y_bands)

# +
#best_catalogue.write("data_tmp/CDFS-SWIRE_cigale_best_extcor_20180129.fits")
n_zphot = 0
for tile in file_list:
        
    file_num = tile.split('_')[4].split('.')[0]
    best_catalogue = Table.read(tile)
    optnir = ((best_catalogue['flag_optnir_det'] == 3) 
          | (best_catalogue['flag_optnir_det'] == 7))
    
    best_catalogue = best_catalogue[optnir]


    best_catalogue = best_catalogue[~np.isnan(best_catalogue["redshift"])]
    n_zphot += len(best_catalogue)
    

    print("Number of sources with z-spec in file {}: {}".format(file_num, len(best_catalogue)))

    # Correction for galactic extinction
    best_catalogue = correct_galactic_extinction(best_catalogue, inplace=True)

    # Convertion to CIGALE format

    best_catalogue = convert_table_for_cigale(best_catalogue, inplace=True, remove_zerofluxes=True)

    remove_unneeded_fluxes(g_bands)
    remove_unneeded_fluxes(u_bands)
    remove_unneeded_fluxes(r_bands)
    remove_unneeded_fluxes(i_bands)
    remove_unneeded_fluxes(z_bands)
    remove_unneeded_fluxes(y_bands)


    best_catalogue.write("data_tmp/tiles/SGP_cigale_optnir_extcor_{}_{}.fits".format(SUFFIX, file_num), overwrite=True)
    
print("Number of zphot optnir objects: {}".format(n_zphot))
# -

# # Catalogue using spectroscopic redshift

# +
n_zspec = 0
for tile in file_list:
        
    file_num = tile.split('_')[4].split('.')[0]
    best_catalogue = Table.read(tile)
    optnir = ((best_catalogue['flag_optnir_det'] == 3) 
          | (best_catalogue['flag_optnir_det'] == 7))
    
    best_catalogue = best_catalogue[optnir]


    best_catalogue = best_catalogue[~np.isnan(best_catalogue["zspec"])]
    best_catalogue.remove_column('redshift')
    best_catalogue["zspec"].name = "redshift"
    n_zspec += len(best_catalogue)

    print("Number of sources with z-spec in file {}: {}".format(file_num, len(best_catalogue)))

    # Correction for galactic extinction
    best_catalogue = correct_galactic_extinction(best_catalogue, inplace=True)

    # Convertion to CIGALE format

    best_catalogue = convert_table_for_cigale(best_catalogue, inplace=True, remove_zerofluxes=True)

    remove_unneeded_fluxes(g_bands)
    remove_unneeded_fluxes(u_bands)
    remove_unneeded_fluxes(r_bands)
    remove_unneeded_fluxes(i_bands)
    remove_unneeded_fluxes(z_bands)
    remove_unneeded_fluxes(y_bands)


    best_catalogue.write("data_tmp/tiles/SGP_cigale_optnir_extcor_zspec_{}_{}.fits".format(SUFFIX, file_num), 
                         overwrite=True)
    
print("Number of zspec optnir objects: {}".format(n_zspec))

StarcoderdataJupyterScriptsDedupFiltered