diff --git a/Metallicity_Stack_Commons/__init__.py b/Metallicity_Stack_Commons/__init__.py index 439b286..d9874d8 100644 --- a/Metallicity_Stack_Commons/__init__.py +++ b/Metallicity_Stack_Commons/__init__.py @@ -5,11 +5,15 @@ import getpass import numpy as np -version = "0.5.0" +version = "0.6.0" lambda0 = [3726.18, 4101.73, 4340.46, 4363.21, 4861.32, 4958.91, 5006.84] line_type = ['Oxy2', 'Balmer', 'Balmer', 'Single', 'Balmer', 'Single', 'Single'] -line_name = ['OII_3727', 'HDELTA', 'HGAMMA', 'OIII_4363', 'HBETA', 'OIII_4958', 'OIII_5007'] +line_name = ['OII_3727', 'HDELTA', 'HGAMMA', 'OIII_4363', 'HBETA', 'OIII_4958', + 'OIII_5007'] + +line_name_short = {"OII": line_name[0], "4363": line_name[3], + "HB": line_name[4], "OIII": line_name[-1]} fitting_lines_dict = {"lambda0": lambda0, "line_type": line_type, "line_name": line_name} diff --git a/Metallicity_Stack_Commons/attenuation.py b/Metallicity_Stack_Commons/analysis/attenuation.py similarity index 98% rename from Metallicity_Stack_Commons/attenuation.py rename to Metallicity_Stack_Commons/analysis/attenuation.py index 9b7dd37..ebb857b 100644 --- a/Metallicity_Stack_Commons/attenuation.py +++ b/Metallicity_Stack_Commons/analysis/attenuation.py @@ -2,7 +2,7 @@ from astropy.table import Table import numpy as np -from . import k_dict +from .. import k_dict HgHb_CaseB = 0.468 # Hg/Hb ratio for zero reddening diff --git a/Metallicity_Stack_Commons/analysis/composite_indv_detect.py b/Metallicity_Stack_Commons/analysis/composite_indv_detect.py index 9e75238..355483c 100644 --- a/Metallicity_Stack_Commons/analysis/composite_indv_detect.py +++ b/Metallicity_Stack_Commons/analysis/composite_indv_detect.py @@ -5,7 +5,7 @@ from astropy.io import ascii as asc from astropy.table import Table -from ..temp_metallicity_calc import metallicity_calculation +from .temp_metallicity_calc import metallicity_calculation from .. import OIII_r from ..column_names import bin_names0, indv_names0, temp_metal_names0 from ..column_names import filename_dict @@ -13,6 +13,11 @@ ID_name = indv_names0[0] bin_ID_name = bin_names0[0] +logR23_name = indv_names0[1] +logO32_name = indv_names0[2] +two_beta_name = indv_names0[5] +three_beta_name = indv_names0[6] + def main(fitspath, dataset, revised=False, det3=True): """ @@ -57,6 +62,13 @@ def main(fitspath, dataset, revised=False, det3=True): bin_id = composite_table['bin_ID'].data bin_temp = composite_table['T_e'].data + # Read in validation table + valid_file = join(fitspath, dataset, filename_dict['bin_valid']) + if not exists(valid_file): + print("ERROR: File not found! "+valid_file) + return + valid_table = asc.read(valid_file) + # Define [indv_em_line_file] indv_em_line_file = join(fitspath, dataset, filename_dict['indv_prop']) if not exists(indv_em_line_file): @@ -78,34 +90,36 @@ def main(fitspath, dataset, revised=False, det3=True): # Populate composite temperature for individual galaxies adopted_temp = np.zeros(len(indv_em_line_table)) bin_id_indv = np.zeros(len(indv_em_line_table)) - for comp_bin, comp_temp in zip(bin_id, bin_temp): + detect_indv = np.zeros(len(indv_em_line_table)) + for comp_bin, comp_temp, detect in zip(bin_id, bin_temp, valid_table['Detection']): bin_idx = np.where(indv_bin_info_table['bin_ID'].data == comp_bin)[0] adopted_temp[bin_idx] = comp_temp bin_id_indv[bin_idx] = comp_bin + detect_indv[bin_idx] = detect + + O2 = indv_em_line_table['OII_3727_Flux_Gaussian'].data # [OII]3726,3728 fluxes + O3 = indv_em_line_table['OIII_5007_Flux_Gaussian'].data # [OIII]5007 fluxes + O3 = O3 * (1+1/OIII_r) # Scale to include OIII4959; Assume 3.1:1 ratio + Hb = indv_em_line_table['HBETA_Flux_Gaussian'].data # H-beta fluxes - O2 = indv_em_line_table['OII_3727_Flux_Gaussian'].data # [OII]3726,3728 fluxes - O3 = indv_em_line_table['OIII_5007_Flux_Gaussian'].data * OIII_r # [OIII]4959,5007 fluxes (Assume 3.1:1 ratio) - Hb = indv_em_line_table['HBETA_Flux_Gaussian'].data # H-beta fluxes + two_beta = O2/Hb + three_beta = O3/Hb + logR23 = np.log10(two_beta + three_beta) + logO32 = np.log10(O3/O2) - if det3: - com_O_log, metal_dict = metallicity_calculation(adopted_temp, O2/Hb, O3/Hb) + if not det3: + metal_dict = metallicity_calculation(adopted_temp, O2/Hb, O3/Hb) else: - det3 = np.where((indv_bin_info_table['Detection'] == 1.0) | (indv_bin_info_table['Detection'] == 0.5))[0] - temp_com_O_log, temp_metal_dict = \ - metallicity_calculation(adopted_temp[det3], O2[det3]/Hb[det3], - O3[det3]/Hb[det3]) - com_O_log = np.zeros(len(indv_em_line_table)) - com_O_log[det3] = temp_com_O_log - - metal_dict = dict() - for key0 in temp_metal_dict.keys(): - metal_dict[key0] = np.zeros(len(indv_em_line_table)) - metal_dict[key0][det3] = temp_metal_dict[key0] + det3_idx = np.where((detect_indv == 1.0) | (detect_indv == 0.5))[0] + metal_dict = \ + metallicity_calculation(adopted_temp, O2/Hb, O3/Hb, det3=det3_idx) # Define [indv_derived_prop_table] to include ID, bin_ID, composite T_e, # and 12+log(O/H) - arr0 = [indv_em_line_table[ID_name], bin_id_indv, adopted_temp, com_O_log] - names0 = [ID_name, bin_ID_name] + temp_metal_names0[:2] + arr0 = [indv_em_line_table[ID_name], bin_id_indv, logR23, logO32, + two_beta, three_beta, adopted_temp] + names0 = [ID_name, bin_ID_name, logR23_name, logO32_name, two_beta_name, + three_beta_name, temp_metal_names0[0]] # Include other metallicities arr0 += list(metal_dict.values()) diff --git a/Metallicity_Stack_Commons/analysis/error_prop.py b/Metallicity_Stack_Commons/analysis/error_prop.py new file mode 100644 index 0000000..69a020f --- /dev/null +++ b/Metallicity_Stack_Commons/analysis/error_prop.py @@ -0,0 +1,161 @@ +from os.path import join, exists + +from chun_codes import random_pdf, compute_onesig_pdf +from . import line_name +from astropy.io import ascii as asc +import numpy as np + +from ..column_names import filename_dict, temp_metal_names0, npz_filename_dict +from .ratios import error_prop_flux_ratios +from .temp_metallicity_calc import temp_calculation, metallicity_calculation + + +def write_npz(path, npz_files, dict_list): + """ + Purpose: + Write numpy files with provided dictionaries + + :param path: str - prefix for filename output + :param npz_files: list - contains npz file names + :param dict_list: list - contains dictionaries for each corresponding npz file + + :return: Write npz files + """ + for file, dict_input in zip(npz_files, dict_list): + npz_outfile = join(path, file) + if exists(npz_outfile): + print("Overwriting : "+npz_outfile) + else: + print("Writing : "+npz_outfile) + np.savez(npz_outfile, **dict_input) + + +def fluxes_derived_prop(path, binned_data=True): + """ + Purpose: + Use measurements and their uncertainties to perform a randomization + approach. The randomization is performed on individual emission lines. + It carries that information to derived flux ratios and then + determines electron temperature and metallicity + + :param path: str of full path + :param binned_data: bool for whether to analysis binned data. Default: True + + :return: + """ + + # Define files to read in for binned data + if binned_data: + flux_file = join(path, filename_dict['bin_fit']) + prop_file = join(path, filename_dict['bin_derived_prop']) + verify_file = join(path, filename_dict['bin_valid']) + + flux_tab0 = asc.read(flux_file) + prop_tab0 = asc.read(prop_file) + + if binned_data: + verify_tab = asc.read(verify_file) + detection = verify_tab['Detection'].data + + # For now we are only considering those with reliable detection and + # excluding those with reliable non-detections (detect = 0.5) + detect_idx = np.where((detection == 1))[0] + + ID = verify_tab['bin_ID'].data + ID_detect = ID[detect_idx] + print(ID_detect) + + flux_tab = flux_tab0[detect_idx] + prop_tab = prop_tab0[detect_idx] + + flux_cols = [str0+'_Flux_Gaussian' for str0 in line_name] + flux_rms_cols = [str0+'_RMS' for str0 in line_name] + + # + # EMISSION-LINE SECTION + # + + # Initialize dictionaries + flux_pdf_dict = dict() + flux_peak = dict() + flux_lowhigh = dict() + + # Randomization for emission-line fluxes + for aa, flux, rms in zip(range(len(flux_cols)), flux_cols, flux_rms_cols): + flux_pdf = random_pdf(flux_tab[flux], flux_tab[rms], seed_i=aa, + n_iter=1000) + err, peak = compute_onesig_pdf(flux_pdf, flux_tab[flux], usepeak=True) + + # Fill in dictionary + flux_pdf_dict[line_name[aa]] = flux_pdf + flux_peak[line_name[aa] + '_peak'] = peak + flux_lowhigh[line_name[aa] + '_lowhigh_error'] = err + + # Update values + flux_tab0[line_name[aa] + '_Flux_Gaussian'][detect_idx] = peak + + # Write revised emission-line fit ASCII table + new_flux_file = join(path, filename_dict['bin_fit_rev']) + if exists(new_flux_file): + print("Overwriting: "+new_flux_file) + else: + print("Writing: "+new_flux_file) + asc.write(flux_tab0, new_flux_file, overwrite=True, format='fixed_width_two_line') + + # Save flux npz files + npz_files = [npz_filename_dict['flux_pdf'], + npz_filename_dict['flux_errors'], + npz_filename_dict['flux_peak']] + dict_list = [flux_pdf_dict, flux_lowhigh, flux_peak] + write_npz(path, npz_files, dict_list) + + # Obtain distributions of line ratios: logR23, logO32, two_beta, three_beta, R + flux_ratios_dict = error_prop_flux_ratios(flux_pdf_dict) + + # + # DERIVED PROPERTIES SECTION + # + + # Initialize dictionaries + derived_prop_pdf_dict = dict() + derived_prop_error_dict = dict() + derived_prop_peak_dict = dict() + + # Calculate temperature distribution + Te_pdf = temp_calculation(flux_ratios_dict['R']) + derived_prop_pdf_dict[temp_metal_names0[0]] = Te_pdf + + # Calculate metallicity distribution + metal_dict = metallicity_calculation(Te_pdf, flux_ratios_dict['two_beta'], + flux_ratios_dict['three_beta']) + derived_prop_pdf_dict.update(metal_dict) + + # Loop for each derived properties (T_e, metallicity, etc.) + for names0 in temp_metal_names0: + arr0 = prop_tab[names0].data + + err_prop, peak_prop = \ + compute_onesig_pdf(derived_prop_pdf_dict[names0], arr0, usepeak=True) + + # Fill in dictionary + derived_prop_error_dict[names0+'_lowhigh_error'] = err_prop + derived_prop_peak_dict[names0+'_peak'] = peak_prop + + # Update values + prop_tab0[names0][detect_idx] = peak_prop + + # Write revised properties ASCII table + new_prop_file = join(path, filename_dict['bin_derived_prop_rev']) + if exists(new_prop_file): + print("Overwriting: "+new_prop_file) + else: + print("Writing: "+new_prop_file) + asc.write(prop_tab0, new_prop_file, overwrite=True, format='fixed_width_two_line') + + # Save derived properties npz files + npz_files = [npz_filename_dict['der_prop_pdf'], + npz_filename_dict['der_prop_errors'], + npz_filename_dict['der_prop_peak']] + dict_list = [derived_prop_pdf_dict, derived_prop_error_dict, + derived_prop_peak_dict] + write_npz(path, npz_files, dict_list) diff --git a/Metallicity_Stack_Commons/fitting.py b/Metallicity_Stack_Commons/analysis/fitting.py similarity index 99% rename from Metallicity_Stack_Commons/fitting.py rename to Metallicity_Stack_Commons/analysis/fitting.py index fcd701f..98e4487 100644 --- a/Metallicity_Stack_Commons/fitting.py +++ b/Metallicity_Stack_Commons/analysis/fitting.py @@ -2,7 +2,7 @@ from astropy.convolution import Box1DKernel, convolve from astropy.io import ascii as asc -from . import scalefact, wavelength_dict +from .. import scalefact, wavelength_dict con1 = wavelength_dict['OII_3729'] / wavelength_dict['OII_3726'] # Ratio of OII doublet line diff --git a/Metallicity_Stack_Commons/analysis/ratios.py b/Metallicity_Stack_Commons/analysis/ratios.py new file mode 100644 index 0000000..227b64e --- /dev/null +++ b/Metallicity_Stack_Commons/analysis/ratios.py @@ -0,0 +1,44 @@ +import numpy as np + +from . import line_name_short, OIII_r +from .temp_metallicity_calc import R_calculation + + +def error_prop_flux_ratios(flux_dict, EBV=None): + """ + Purpose: + Primary code to determine a variety of line ratios based on a dictionary + containing emission-line fluxes + + :param flux_dict: dictionary containing line ratios + :param EBV: array of E(B-V). Same dimensions as contents of flux_dict + + :return flux_ratios_dict: dictionary containing flux ratios + """ + + # Retrieve emission line fluxes + OII = flux_dict[line_name_short['OII']] + OIII = flux_dict[line_name_short['OIII']] + Hb = flux_dict[line_name_short['HB']] + OIII4363 = flux_dict[line_name_short['4363']] + + # Define flux ratios + two_beta = OII/Hb + three_beta = (1+1/OIII_r) * OIII/Hb + logR23 = np.log10(two_beta + three_beta) + logO32 = np.log10((1+1/OIII_r) * OIII/OII) + + # Define dictionary of flux ratios + flux_ratios_dict = dict() + flux_ratios_dict['two_beta'] = two_beta + flux_ratios_dict['three_beta'] = three_beta + flux_ratios_dict['logR23'] = logR23 + flux_ratios_dict['logO32'] = logO32 + + if EBV is None: + print("Not applying dust attenuation correction") + EBV = np.zeros(OIII4363.shape) + + flux_ratios_dict['R'] = R_calculation(OIII4363, OIII, EBV) + + return flux_ratios_dict diff --git a/Metallicity_Stack_Commons/temp_metallicity_calc.py b/Metallicity_Stack_Commons/analysis/temp_metallicity_calc.py similarity index 58% rename from Metallicity_Stack_Commons/temp_metallicity_calc.py rename to Metallicity_Stack_Commons/analysis/temp_metallicity_calc.py index 5c1af77..b4312a1 100644 --- a/Metallicity_Stack_Commons/temp_metallicity_calc.py +++ b/Metallicity_Stack_Commons/analysis/temp_metallicity_calc.py @@ -2,7 +2,7 @@ from . import k_dict -from .column_names import temp_metal_names0, remove_from_list +from ..column_names import temp_metal_names0, remove_from_list k_4363 = k_dict['OIII_4363'] k_5007 = k_dict['OIII_5007'] @@ -48,7 +48,7 @@ def temp_calculation(R): return T_e -def metallicity_calculation(T_e, TWO_BETA, THREE_BETA): +def metallicity_calculation(T_e, TWO_BETA, THREE_BETA, det3=None): """ Determines 12+log(O/H) from electron temperature and [OII]/Hb and [OIII]/Hb flux ratio @@ -56,27 +56,38 @@ def metallicity_calculation(T_e, TWO_BETA, THREE_BETA): :param TWO_BETA: numpy array of [OII]/Hb flux ratio :param THREE_BETA: numpy array of [OIII]/Hb flux ratio - :return com_O_log: numpy array of 12+log(O/H) - :return metal_dict: dictionary containing O+/H, O++/H, log(O+/H), log(O++/H) + :return metal_dict: dictionary containing 12+log(O/H), O+/H, O++/H, log(O+/H), log(O++/H) """ - t_3 = T_e * 1e-4 - t_2 = 0.7 * t_3 + 0.17 - x2 = 1e-4 * 1e3 * t_2 ** (-0.5) + n_sample = len(T_e) + t_3 = np.zeros(n_sample) + t_2 = np.zeros(n_sample) + x2 = np.zeros(n_sample) + + if det3 is None: + det3 = np.arange(n_sample) + + t_3[det3] = T_e[det3] * 1e-4 + t_2[det3] = 0.7 * t_3[det3] + 0.17 + x2[det3] = 1e-4 * 1e3 * t_2[det3] ** (-0.5) + + O_s_ion_log = np.zeros(n_sample) + O_d_ion_log = np.zeros(n_sample) # Equations from Izotov et al. (2006) - O_s_ion_log = np.log10(TWO_BETA) + 5.961 + 1.676 / t_2 - 0.4 * np.log10(t_2) \ - - 0.034 * t_2 + np.log10(1 + 1.35 * x2) - 12 - O_d_ion_log = np.log10(THREE_BETA) + 6.200 + 1.251 / t_3 \ - - 0.55 * np.log10(t_3) - 0.014 * (t_3) - 12 + O_s_ion_log[det3] = np.log10(TWO_BETA[det3]) + 5.961 + 1.676 / t_2[det3] \ + - 0.4 * np.log10(t_2[det3]) - 0.034 * t_2[det3] + \ + np.log10(1 + 1.35 * x2[det3]) - 12 + O_d_ion_log[det3] = np.log10(THREE_BETA[det3]) + 6.200 + 1.251 / t_3[det3] \ + - 0.55 * np.log10(t_3[det3]) - 0.014 * (t_3[det3]) - 12 O_s_ion = 10 ** O_s_ion_log O_d_ion = 10 ** O_d_ion_log com_O = O_s_ion + O_d_ion com_O_log = np.log10(com_O) + 12 - key_dict = remove_from_list(temp_metal_names0, ['T_e', '12+log(O/H)']) - key_values = [O_s_ion_log, O_d_ion_log, O_s_ion, O_d_ion] # Order matters here + key_dict = remove_from_list(temp_metal_names0, ['T_e']) + key_values = [com_O_log, O_s_ion_log, O_d_ion_log, O_s_ion, O_d_ion] # Order matters here metal_dict = dict(zip(key_dict, key_values)) - return com_O_log, metal_dict + return metal_dict diff --git a/Metallicity_Stack_Commons/column_names.py b/Metallicity_Stack_Commons/column_names.py index 62bb685..f191615 100644 --- a/Metallicity_Stack_Commons/column_names.py +++ b/Metallicity_Stack_Commons/column_names.py @@ -82,6 +82,15 @@ def line_fit_suffix_add(line_name0, line_type0): filename_dict['indv_bin_info'] = 'individual_bin_info.tbl' filename_dict['indv_derived_prop'] = 'individual_derived_properties.tbl' +# numpy files +npz_filename_dict = dict() +npz_filename_dict['flux_pdf'] = 'flux_pdf.npz' +npz_filename_dict['flux_errors'] = 'flux_errors.npz' +npz_filename_dict['flux_peak'] = 'flux_peak.npz' +npz_filename_dict['der_prop_pdf'] = 'derived_properties_pdf.npz' +npz_filename_dict['der_prop_errors'] = 'derived_properties_errors.npz' +npz_filename_dict['der_prop_peak'] = 'derived_properties_peak.npz' + def merge_column_names(*args): """ diff --git a/Metallicity_Stack_Commons/plotting/balmer.py b/Metallicity_Stack_Commons/plotting/balmer.py index e042313..7addc0e 100644 --- a/Metallicity_Stack_Commons/plotting/balmer.py +++ b/Metallicity_Stack_Commons/plotting/balmer.py @@ -14,7 +14,7 @@ from astropy.io import ascii as asc from matplotlib.backends.backend_pdf import PdfPages -from ..fitting import gauss, double_gauss +from ..analysis.fitting import gauss, double_gauss from .. import scalefact, wavelength_dict diff --git a/setup.py b/setup.py index 14d774f..c0e6c3b 100644 --- a/setup.py +++ b/setup.py @@ -5,7 +5,7 @@ setup( name='Metallicity_Stack_Commons', - version='0.5.0', + version='0.6.0', packages=find_packages('Metallicity_Stack_Commons'), url='https://github.com/astrochun/Metallicity_Stack_Commons', license='MIT License',