Last updated: Bryce Kalmbach, December 2018
Danny Goldstein's latest strongly lensed SNe catalog has all the properties needed to replace galaxies from DC2 with the sprinkler except for SEDs and CatSim magNorm values for the lens galaxies. This notebook matches the lensed SNe lens galaxies to a DC2 galaxy then uses sims_GCRCatSimInterface and the given 'i-band' magnitude from the lens catalog to calculate an appropriate magNorm value and SED.
You will need gcrCatalogs, GCR, pandas, numpy and matplotlib
In [26]:
import numpy as np
import pandas as pd
import matplotlib as mpl
mpl.rcParams['text.usetex'] = False
from matplotlib import pyplot as plt
import os
from lsst.sims.photUtils import Sed, Bandpass, BandpassDict, getImsimFluxNorm
%matplotlib inline
Set the gcr-catalog version you need below and all names of outputs will be consistent.
In [27]:
catalog_version = 'cosmoDC2_v1.1.4'
In [28]:
dc2_df_system = pd.read_hdf('/global/cscratch1/sd/brycek/glsne_dc2_v2.h5', key='system')
dc2_df_image = pd.read_hdf('/global/cscratch1/sd/brycek/glsne_dc2_v2.h5', key='image')
In [29]:
import GCRCatalogs
import pandas as pd
from GCR import GCRQuery
In [30]:
# _small version of cosmoDC2 catalog will work for this
catalog = GCRCatalogs.load_catalog('%s_small' % catalog_version)
In [31]:
# Example lens catalog row.
dc2_df_system.iloc[0]
Out[31]:
In [8]:
%%time
z_cat_min = np.power(10, np.log10(np.min(dc2_df_system['zl'])) - .1)
z_cat_max = np.power(10, np.log10(np.max(dc2_df_system['zl'])) + .1)
i_gal_dim = np.max(-2.5*np.log10(np.power(10, -0.4*dc2_df_system['lensgal_mi'])*.9))
i_gal_bright = np.min(-2.5*np.log10(np.power(10, -0.4*dc2_df_system['lensgal_mi'])*1.1))
print(z_cat_min, z_cat_max, i_gal_dim, i_gal_bright)
data = catalog.get_quantities(['galaxy_id', 'redshift_true', 'ellipticity_true', 'mag_true_i_lsst',
'size_minor_bulge_true', 'size_bulge_true'],
filters=['mag_true_i_lsst > %f' % i_gal_bright,
'mag_true_i_lsst < %f' % i_gal_dim,
'redshift_true > %f' % z_cat_min, 'redshift_true < %f' % z_cat_max,
'stellar_mass_bulge/stellar_mass > 0.99'])
data['glsne_ellipticity'] = (1-(data['size_minor_bulge_true']/data['size_bulge_true']))
data_df = pd.DataFrame(data)
print(data_df.head(10))
In [9]:
data_df.to_csv('sn_matching_checkpoint_1.csv', index=False)
In [10]:
#data_df = pd.read_csv('sn_matching_checkpoint_1.csv')
In [11]:
%%time
gcr_glsn_match = []
err = 0
np.random.seed(10)
i = 0
row_num = -1
keep_rows = []
num_match_total = []
for zsrc, i_mag_star in zip(dc2_df_system['zl'].values,
dc2_df_system['lensgal_mi'].values):
row_num += 1
#print(zsrc, m_star, ellip)
if row_num % 10000 == 0:
print(row_num)
z_min, z_max = np.power(10, np.log10(zsrc) - .1), np.power(10, np.log10(zsrc) + .1)
i_mag_dim, i_mag_bright = -2.5*np.log10(np.power(10, -0.4*i_mag_star)*.9),\
-2.5*np.log10(np.power(10, -0.4*i_mag_star)*1.1)
#ellip_min, ellip_max = ellip*.9, ellip*1.1
#r_major = r_eff / np.sqrt(1 - ellip)
#r_major_max, r_major_min = r_major*1.1, r_major*.9
#print(z_min, z_max, i_mag_min, i_mag_max, ellip_min, ellip_max)
data_subset = data_df.query('redshift_true > %f and redshift_true < %f ' % (z_min, z_max)
+ 'and mag_true_i_lsst > %f and mag_true_i_lsst < %f ' % (i_mag_bright, i_mag_dim)
#+ 'and glsne_ellipticity > %f and glsne_ellipticity < %f' % (ellip_min, ellip_max)
#+ 'and size_bulge_true > %f and size_bulge_true < %f' % (r_major_min, r_major_max)
)
num_matches = len(data_subset['redshift_true'])
num_match_total.append(num_matches)
if num_matches == 0:
err += 1
continue
elif num_matches == 1:
gcr_data = [data_subset['redshift_true'].values[0],
data_subset['galaxy_id'].values[0]]
gcr_glsn_match.append(gcr_data)
keep_rows.append(row_num)
elif num_matches > 1:
use_idx = np.random.choice(num_matches)
gcr_data = [data_subset['redshift_true'].values[use_idx],
data_subset['galaxy_id'].values[use_idx]]
gcr_glsn_match.append(gcr_data)
keep_rows.append(row_num)
print("Total Match Failures: ", err, " Percentage Match Failures: ", np.float(err)/len(dc2_df_system))
In [12]:
gcr_glsn_match = np.array(gcr_glsn_match)
keep_rows = np.array(keep_rows)
In [13]:
n, bins, _ = plt.hist(dc2_df_system['zl'], alpha=0.4, bins=20)
plt.hist(gcr_glsn_match[:,0], alpha=0.4, bins=bins)
Out[13]:
In [14]:
np.savetxt('gcr_glsn_match.txt', gcr_glsn_match)
np.savetxt('keep_rows_sn.txt', keep_rows)
It seems we are missing mostly very low redshift systems. This is most likely a result of our percentage bounds on redshift. However, we still have nearly 40000 systems in this low redshift range and will have enough to match.
sims_GCRCatSimInterfaceOur goal is to add SEDs and magNorm values into the lensing catalogs for the lens galaxies. Here we get the top hat filters out of cosmoDC2 and use the code in sims_GCRCatSimInterface to match these values to a CATSIM SED file in the same way the galaxies are matched for Instance Catalog production in DC2. We also use the code to calculate the magnitude normalization for PhoSim.
In [7]:
import sys
sys.path.append('/global/homes/b/brycek/DC2/sims_GCRCatSimInterface/workspace/sed_cache/')
In [8]:
from SedFitter import sed_from_galacticus_mags
In [9]:
H0 = catalog.cosmology.H0.value
Om0 = catalog.cosmology.Om0
In [10]:
sed_label = []
sed_min_wave = []
sed_wave_width = []
for quant_label in sorted(catalog.list_all_quantities()):
if (quant_label.startswith('sed') and quant_label.endswith('bulge')):
sed_label.append(quant_label)
label_split = quant_label.split('_')
sed_min_wave.append(int(label_split[1])/10)
sed_wave_width.append(int(label_split[2])/10)
bin_order = np.argsort(sed_min_wave)
sed_label = np.array(sed_label)[bin_order]
sed_min_wave = np.array(sed_min_wave)[bin_order]
sed_wave_width = np.array(sed_wave_width)[bin_order]
Check to see that our bins are now in order when we call them.
In [11]:
for i in zip(sed_label, sed_min_wave, sed_wave_width):
print(i)
In [20]:
del(data)
del(data_df)
In [21]:
%%time
columns = ['galaxy_id', 'redshift_true', 'mag_u_lsst', 'mag_g_lsst', 'mag_r_lsst',
'mag_i_lsst', 'mag_z_lsst', 'mag_y_lsst']
for sed_bin in sed_label:
columns.append(sed_bin)
data = catalog.get_quantities(columns,
filters=['mag_true_i_lsst > %f' % i_gal_bright,
'mag_true_i_lsst < %f' % i_gal_dim,
'redshift_true > %f' % z_cat_min, 'redshift_true < %f' % z_cat_max,
'stellar_mass_bulge/stellar_mass > 0.99'])
data_df = pd.DataFrame(data)
In [22]:
data_df.to_csv('sn_matching_checkpoint_2.csv', index=False)
In [32]:
# _small version of cosmoDC2 catalog will work for this
# catalog = GCRCatalogs.load_catalog('%s_small' % catalog_version)
data_df = pd.read_csv('sn_matching_checkpoint_2.csv')
gcr_glsn_match = np.genfromtxt('gcr_glsn_match.txt')
keep_rows = np.genfromtxt('keep_rows_sn.txt')
In [33]:
keep_rows = np.array(keep_rows, dtype=int)
In [24]:
%%time
sed_name_list = []
magNorm_list = []
i = 0
lsst_mags = []
mag_30_list = []
redshift_list = []
for gal_id, gal_z in zip(gcr_glsn_match[:, 1],
gcr_glsn_match[:, 0]):
if i % 10000 == 0:
print(i)
i+=1
#print(int(gal_id))
data_subset = data_df.query(str('galaxy_id == %i' % int(gal_id))) ## Galaxy Ids are not unique in cosmoDC2_v0.1
#print(data_subset)
mag_array = []
lsst_mag_array = [data_subset['mag_%s_lsst' % band_name].values[0] for band_name in ['u', 'g', 'r', 'i', 'z', 'y']]
for sed_bin in sed_label:
mag_array.append(-2.5*np.log10(data_subset[sed_bin].values[0]))
mag_array = np.array(mag_array)
lsst_mag_array = np.array(lsst_mag_array)
lsst_mags.append(lsst_mag_array)
mag_30_list.append(mag_array)
redshift_list.append(gal_z)
mag_30_list = np.array(mag_30_list).T
lsst_mags = np.array(lsst_mags).T
redshift_list = np.array(redshift_list)
In [25]:
mag_30_list[:,0], lsst_mags[:,0]
Out[25]:
In [26]:
print(np.shape(mag_30_list), np.shape(lsst_mags), np.shape(redshift_list))
In [27]:
sed_name, magNorm, av, rv = sed_from_galacticus_mags(mag_30_list, redshift_list,
H0, Om0, sed_min_wave, sed_wave_width, lsst_mags)
sed_name_list = sed_name
magNorm_list = magNorm
print(len(sed_name_list), len(keep_rows))
In [30]:
sed_name_array = np.array(sed_name_list)
magNorm_array = np.array(magNorm_list)
av_array = np.array(av)
rv_array = np.array(rv)
In [31]:
np.savetxt('sed_name_list_sn.txt', sed_name_array, fmt='%s')
np.savetxt('magNorm_list_sn.txt', magNorm_array)
np.savetxt('av_list_sn.txt', av_array)
np.savetxt('rv_list_sn.txt', rv_array)
In [34]:
sed_name_array = np.genfromtxt('sed_name_list_sn.txt', dtype=np.str)
magNorm_array = np.genfromtxt('magNorm_list_sn.txt')
av_array = np.genfromtxt('av_list_sn.txt')
rv_array = np.genfromtxt('rv_list_sn.txt')
redshift_list = gcr_glsn_match[:,0]
In [35]:
test_bandpassDict = BandpassDict.loadTotalBandpassesFromFiles()
In [36]:
imsimband = Bandpass()
imsimband.imsimBandpass()
In [38]:
mag_norm_glsne = []
test_sed = Sed()
test_sed.readSED_flambda(os.path.join(str(os.environ['SIMS_SED_LIBRARY_DIR']), sed_name_array[0]))
a, b = test_sed.setupCCM_ab()
# We need to adjust the magNorms of the galaxies so that,
# in the i-band, they match the OM10 APMAG_I.
# We do this be calculating the i-magnitudes of the galaxies as they
# will be simulated, finding the difference between that magnitude
# and APMAG_I, and adding that difference to *all* of the magNorms
# of the galaxy (this way, the cosmoDC2 validated colors of the
# galaxies are preserved)
for i, idx in list(enumerate(keep_rows)):
if i % 10000 == 0:
print(i, idx)
test_sed = Sed()
test_sed.readSED_flambda(os.path.join(str(os.environ['SIMS_SED_LIBRARY_DIR']), sed_name_array[i]))
fnorm = getImsimFluxNorm(test_sed, magNorm_array[3,i])
test_sed.multiplyFluxNorm(fnorm)
test_sed.addDust(a, b, A_v=av_array[i], R_v=rv_array[i])
test_sed.redshiftSED(redshift_list[i], dimming=True)
i_mag = test_sed.calcMag(test_bandpassDict['i'])
d_mag = dc2_df_system['lensgal_mi'].iloc[idx]-i_mag
mag_norm_glsne.append(magNorm_array[:,i] + d_mag)
# # Double check colors are maintained
# test_sed = Sed()
# test_sed.readSED_flambda(os.path.join(str(os.environ['SIMS_SED_LIBRARY_DIR']), sed_name_array[i]))
# print(magNorm_array[3,i])
# fn_1 = test_sed.calcFluxNorm(magNorm_array[3,i], imsimband)
# test_sed.multiplyFluxNorm(fn_1)
# i_orig = test_sed.calcMag(test_bandpassDict['i'])
# fn_2 = test_sed.calcFluxNorm(magNorm_array[2,i], imsimband)
# test_sed.multiplyFluxNorm(fn_2)
# r_orig = test_sed.calcMag(test_bandpassDict['r'])
# print(r_orig, i_orig, r_orig-i_orig)
# fn_3 = test_sed.calcFluxNorm(mag_norm_glsne[-1][2], imsimband)
# test_sed.multiplyFluxNorm(fn_3)
# r_new = test_sed.calcMag(test_bandpassDict['r'])
# print(dc2_df_system['lensgal_mi'].iloc[idx], r_new, r_new - dc2_df_system['lensgal_mi'].iloc[idx])
# print(mag_norm_glsne[-1])
In [39]:
mag_norm_glsne = np.array(mag_norm_glsne)
In [40]:
sed_metals = []
sed_ages = []
sed_metals_dict = {'0005Z':'.005', '002Z':'.02', '02Z':'.2', '04Z':'.4', '1Z':'1.0', '25Z':'2.5'}
for sed_template in sed_name_array:
sed_info = sed_template.split('/')[1].split('.')
sed_age_info = sed_info[1].split('E')
sed_ages.append(np.power(10, int(sed_age_info[1]))*int(sed_age_info[0]))
sed_metals.append(sed_metals_dict[sed_info[2]])
In [41]:
fig = plt.figure(figsize=(12, 6))
mpl.rcParams['text.usetex'] = False
fig.add_subplot(1,2,1)
plt.hist(np.log10(sed_ages))
plt.xlabel('Log10(Galaxy Age (years))')
plt.ylabel('Counts')
fig.add_subplot(1,2,2)
names, counts = np.unique(sed_metals, return_counts=True)
x = np.arange(len(names), dtype=int)
plt.bar(x, counts)
plt.xticks(x, names)
plt.xlabel('Metallicity (Z_sun)')
plt.ylabel('Counts')
plt.tight_layout()
plt.subplots_adjust(top=0.9)
plt.suptitle('Age and Metallicity of matched SED templates')
Out[41]:
We want our lens galaxies to be large, red ellipticals and it looks like we are getting older, redder galaxies which looks like we are succeeding.
We will take all the columns currently in the Twinkles GLSNe catalog and add in our new SED filenames and SED magnitude normalizations for PhoSim.
In [42]:
# Double check the distribution of t0
mpl.rcParams['text.usetex'] = False
plt.hist(dc2_df_system['t0'])
plt.xlabel('Start Date (MJD)')
plt.ylabel('Count')
Out[42]:
In the final catalog we produce we only want to include the systems we were able to match to a gcr galaxy.
In [43]:
results_dict = {}
print(len(keep_rows), len(gcr_glsn_match[:, 0]))
for keep_idx in range(len(keep_rows)):
results_dict[str(dc2_df_system['sysno'].iloc[keep_rows[keep_idx]])] = {'z':gcr_glsn_match[:, 0][keep_idx],
'sed_name':sed_name_array[keep_idx],
'magNorm':mag_norm_glsne[keep_idx],
'lens_av':av_array[keep_idx],
'lens_rv':rv_array[keep_idx]}
In [44]:
keep_systems = dc2_df_system['sysno'].iloc[keep_rows].values
In [45]:
final_df_z = []
final_df_lens_sed = []
final_df_magNorm_u = []
final_df_magNorm_g = []
final_df_magNorm_r = []
final_df_magNorm_i = []
final_df_magNorm_z = []
final_df_magNorm_y = []
final_df_lens_av = []
final_df_lens_rv = []
keep_in_df = []
for idx, twinkles_sys in enumerate(dc2_df_system['sysno']):
if twinkles_sys in keep_systems:
keep_in_df.append(idx)
final_df_z.append(results_dict[str(twinkles_sys)]['z'])
final_df_lens_sed.append(results_dict[str(twinkles_sys)]['sed_name'])
final_df_magNorm_u.append(results_dict[str(twinkles_sys)]['magNorm'][0])
final_df_magNorm_g.append(results_dict[str(twinkles_sys)]['magNorm'][1])
final_df_magNorm_r.append(results_dict[str(twinkles_sys)]['magNorm'][2])
final_df_magNorm_i.append(results_dict[str(twinkles_sys)]['magNorm'][3])
final_df_magNorm_z.append(results_dict[str(twinkles_sys)]['magNorm'][4])
final_df_magNorm_y.append(results_dict[str(twinkles_sys)]['magNorm'][5])
final_df_lens_av.append(results_dict[str(twinkles_sys)]['lens_av'])
final_df_lens_rv.append(results_dict[str(twinkles_sys)]['lens_rv'])
In [46]:
final_df = dc2_df_system.iloc[keep_in_df]
final_df = final_df.reset_index(drop=True)
In [47]:
final_df['lensgal_magnorm_u'] = final_df_magNorm_u
final_df['lensgal_magnorm_g'] = final_df_magNorm_g
final_df['lensgal_magnorm_r'] = final_df_magNorm_r
final_df['lensgal_magnorm_i'] = final_df_magNorm_i
final_df['lensgal_magnorm_z'] = final_df_magNorm_z
final_df['lensgal_magnorm_y'] = final_df_magNorm_y
final_df['lensgal_sed'] = final_df_lens_sed
final_df['lens_av'] = final_df_lens_av
final_df['lens_rv'] = final_df_lens_rv
In [48]:
final_df.head()
Out[48]:
Great! Looks like everything is in our new catalog and we can now save it to file.
In [49]:
final_df.to_hdf('/global/cscratch1/sd/brycek/glsne_%s.h5' % catalog_version, key='system', format='table')
In [50]:
dc2_df_image.to_hdf('/global/cscratch1/sd/brycek/glsne_%s.h5' % catalog_version, mode='a', key='image', format='table')
In [ ]: