#!/usr/bin/env python3 import gc import os import sys import gzip import glob import time import shutil import itertools import numpy as np import configparser from datetime import datetime from os import makedirs, remove from scipy.spatial import KDTree from copy import deepcopy as copy from os.path import basename, isfile from scipy.ndimage import median_filter from contextlib import redirect_stdout, redirect_stderr from pyFIT3D.common.io import ReadArguments, print_done, print_block_init from pyFIT3D.common.io import get_data_from_fits, get_wave_from_header, print_time from pyFIT3D.common.io import trim_waves, output_spectra, remove_isfile, write_img_header from pyFIT3D.common.tools import imarith, get_slice, clean_nan, indices_spec from pyFIT3D.common.tools import smooth_spec_clip_cube, map_auto_ssp_rnd_seg from pyFIT3D.common.tools import flux_elines_cube_EW, pack_NAME, sum_mass_age from pyFIT3D.common.tools import spec_extract_cube_mean, spec_extract_cube_error from pyFIT3D.common.tools import radial_sum_cube_e, cont_seg_all_SN, rss_seg2cube # from pyFIT3D.common.tools import get_index_output_auto_ssp_elines_rnd_rss_outFIT3D from pyFIT3D.common.tools import read_img_header, redshift_config from pyFIT3D.common.tools import get_index, get_cosmo, pack_results_name from pyFIT3D.common.tools import index_seg_cube, clean_Ha_map, med2df, array_to_fits from pyFIT3D.common.tools import vel_eline, SN_map_seg, get_SN_cube, get_SN_rss from pyFIT3D.common.auto_ssp_tools import ConfigAutoSSP from pyFIT3D.common.auto_ssp_tools import auto_ssp_elines_rnd_rss_main from pyFIT3D.common.auto_ssp_tools import auto_ssp_elines_single_main, dump_rss_output from pyFIT3D.common.gas_tools import output_emission_lines_spectra, create_emission_lines_parameters from pyFIT3D.common.gas_tools import kin_cube_elines_main, output_emission_lines_parameters, append_emission_lines_parameters from pyFIT3D.common.constants import __selected_R_V__, __selected_extlaw__, __n_Monte_Carlo__ from pyFIT3D.common.constants import __c__, __Ha_central_wl__, _MODELS_ELINE_PAR, _SSP_OUTPUT_INDEX from pyFIT3D.common.constants import __FWHM_to_sigma__, __sigma_to_FWHM__, __indices__, __solar_metallicity__ from pyFIT3D.common.constants import __Hubble_constant__, __Omega_matter__, __Omega_Lambda__, __solar_luminosity__ from pyFIT3D.modelling.stellar import SSPModels, StPopSynt class ReadArgumentsLocal(ReadArguments): """ Argument parser for this script """ __script_name__ = basename(sys.argv[0]) __mandatory__ = ['name', 'config_file'] __optional__ = ['redshift', 'x0', 'y0', 'only_binned_spaxels'] __optional__ = ['redshift', 'x0', 'y0'] #, 'only_binned_spaxels'] __arg_names__ = __mandatory__ + __optional__ __N_tot_args__ = len(__arg_names__) __conv_func_mandatory__ = {'name': str, 'config_file': str} # __conv_func_optional__ = {'only_binned_spaxels': bool} __conv_func__ = __conv_func_mandatory__.copy() # __conv_func__.update(__conv_func_optional__) usage_msg_tmp1 = 'USE: {}'.format(__script_name__) usage_msg_tmp2 = ' NAME CONFIG_FILE' usage_msg_tmp3 = ' [REDSHIFT] [X0] [Y0] [ONLY_BINNED_SPAXELS=0/1]' usage_msg_tmp4 = '\nCONFIG_FILE is mandatory but defaults to ana_single.ini' # usage_msg_tmp5 = '\nONLY_BINNED_SPAXELS=1 will reanalyze only spaxels in voxels with more than 1 binned spaxel.' # usage_msg_tmp6 = '\n\tDefaults to 0, i.e., reanalyze all spaxels associated to voxels.' __usage_msg__ = usage_msg_tmp1 + usage_msg_tmp2 + usage_msg_tmp3 + usage_msg_tmp4 # __usage_msg__ += usage_msg_tmp5 + usage_msg_tmp6 def __init__(self, args_list=None, verbose=False): ReadArguments.__init__(self, args_list, verbose=verbose) self._parse() def _parse(self): if self.config_file is None: self.config_file = 'ana_single.ini' if not isfile(self.config_file): print(f'{self.__script_name__}: {self.config_file}: not found') sys.exit() # EL: TODO def get_redshift_from_file(filename): z = None if filename is None: return None # if isfile(filename): # Read the file FILENAME and get redshift. return z # XXX: EADL # It should be called mask_regions_brighter_than_center() def mask_stars_in_FoV(name, V__yx, msk__yx, x0=None, y0=None, mask_val=0): time_ini_loc = print_block_init('Mask stars in FoV...') ny, nx = V__yx.shape ix0 = int(nx*.5) if x0 is None else x0 iy0 = int(ny*.5) if y0 is None else y0 V_msk__yx = np.ones_like(msk__yx) V_msk__yx[msk__yx == 2] = 0 # look for stars in the field ix0, iy0 = int(nx*.5), int(ny*.5) val0 = V__yx[iy0, ix0] max_dist = 0.15*np.sqrt(nx**2 + ny**2) for ixy in itertools.product(range(nx),range(ny)): ix, iy = ixy val = V__yx[iy, ix] dist = np.sqrt((ix - ix0)**2 + (iy - iy0)**2) if (val > val0) & (dist > max_dist): nx0 = ix - 3 nx1 = ix + 3 ny0 = iy - 3 ny1 = iy + 3 # img border control nx0 = 0 if nx0 < 0 else nx0 ny0 = 0 if ny0 < 0 else ny0 nx1 = nx - 1 if nx1 > (nx - 1) else nx1 ny1 = ny - 1 if ny1 > (ny - 1) else ny1 # mask region for iixy in itertools.product(range(nx0, nx1),range(ny0, ny1)): iix, iiy = iixy V_msk__yx[iiy, iix] = mask_val print_done(time_ini=time_ini_loc) return V_msk__yx def collapse_badpixels_mask(mask__wyx, threshold_fraction=1, mask_value=1): """ Performs the bad pixels mask collapse, i.e., generates a 2D map from a 3D (wavelengths, y, x) masking spaxels with a fraction of masked pixels equal or above `threshold`. Parameters ---------- mask__wyx : array like Bad pixels mask cube (3 dimensions, NWAVE, NY, NX). threshold_fraction : float Sets the threshold to be considered a bad spaxel, i.e. spaxels with a fraction of masked pixels equal or above this threshold will be masked. threshold = NWAVE * threshold_fraction. Returns ------- array like 2D bad spaxels map. """ nw, ny, nx = mask__wyx.shape threshold = nw*threshold_fraction badspaxels__yx = np.zeros((ny, nx), dtype=mask__wyx.dtype) badpixels_total__yx = mask__wyx.sum(axis=0) badspaxels__yx[badpixels_total__yx/nw >= threshold] == mask_value return badspaxels__yx def create_SSP_from_spectra(ssp_library, coeffs, redshift=0, sigma=0, AV=0, sigma_inst=0, RV=3.1, extlaw='CCM', output_filename=None): nm = coeffs.shape[-1] # CREATING MODELS flux_models__mw = np.array( [ssp_library.get_model_from_coeffs( coeffs=coeffs[:, i], wavelength=ssp_library.wavelength, redshift=redshift, sigma=sigma, sigma_inst=sigma_inst, AV=AV, R_V=RV, extlaw=extlaw) for i in range(nm)] ) # CREATING HEADER FROM DICT hdr = {} for i in range(nm): l_age_LW, l_met_LW, _, _ = ssp_library.get_tZ_from_coeffs(coeffs[:, i], mean_log=True) ML = np.dot(coeffs[:, i], ssp_library.mass_to_light) age_LW = 10**l_age_LW met_LW = 10**l_met_LW _strage = f'{age_LW:0.4f}Gyr' _strZ = f'{met_LW:0.7f}'.replace('0.', 'z') hdr[f'NAME{i}'] = f'spec_ssp_{_strage}_{_strZ}.dat' hdr[f'NORM{i}'] = 1/ML hdr['WAVENORM'] = ssp_library.normalization_wavelength hdr['CDELT1'] = ssp_library._header.get('CDELT1') hdr['CRPIX1'] = ssp_library._header.get('CRPIX1') hdr['CRVAL1'] = ssp_library._header.get('CRVAL1') # SAVE FILE USING pyPipe3D TAG! if output_filename is not None: array_to_fits(output_filename, flux_models__mw, header=hdr, overwrite=True) return flux_models__mw def create_deseg_packs(name, org_SFH_pack, org_SSP_pack, org_ELINES_pack, org_output_prefix, org_cont_seg, org_scale_seg, prefix): new_SFH_pack = f'{prefix}.{name}.{basename(org_SFH_pack)}' new_SSP_pack = f'{prefix}.{name}.{basename(org_SSP_pack)}' new_ELINES_pack = f'{prefix}.{name}.{basename(org_ELINES_pack)}' new_output_prefix = org_output_prefix.replace('CS', prefix) f_out = open(new_SFH_pack, 'w') with open(org_SFH_pack, 'r') as f: for l in f.readlines(): if org_output_prefix in l: l = l.replace(org_output_prefix, new_output_prefix) f_out.write(l) f_out.close() f_out = open(new_ELINES_pack, 'w') with open(org_ELINES_pack, 'r') as f: for l in f.readlines(): if org_output_prefix in l: l = l.replace('NAME', f'{prefix}.NAME') f_out.write(l) f_out.close() f_out = open(new_SSP_pack, 'w') with open(org_SSP_pack, 'r') as f: for l in f.readlines(): if org_output_prefix in l: l = l.replace(org_output_prefix, new_output_prefix) elif org_cont_seg in l: l = l.replace(org_cont_seg, f'{prefix}.{org_cont_seg}') elif (org_scale_seg in l): l = l.replace(org_scale_seg, f'{prefix}.{org_scale_seg}') f_out.write(l) f_out.close() return new_SFH_pack, new_SSP_pack, new_ELINES_pack def get_index_cube(wave__w, flux_ssp__wyx, res__wyx, redshift__yx, n_sim, plot=0, wl_half_range=200): nw, ny, nx = flux_ssp__wyx.shape names__i = list(__indices__.keys()) indices = {} for name in names__i: indices[name] = { 'EW': np.zeros((ny, nx), dtype='float'), 'sigma_EW': np.zeros((ny, nx), dtype='float') } indices['SN'] = np.zeros((ny, nx), dtype='float') indices['e_SN'] = np.zeros((ny, nx), dtype='float') k_list = [x for x in indices.keys() if 'SN' not in x] nI = len(k_list) output_cube__Iyx = np.zeros((nI*2 + 2, ny, nx)) for ixy in itertools.product(range(nx),range(ny)): ix, iy = ixy redshift = redshift__yx[iy, ix] res__w = res__wyx[:, iy, ix] flux_ssp__w = flux_ssp__wyx[:, iy, ix] + res__w EW__ki, med_flux, std_flux = indices_spec(wave__w, flux_ssp__w, res__w, redshift, n_sim, plot=plot, wl_half_range=200) for ii, name in enumerate(names__i): # in range(n_ind): # name = names__i[ii] indices[name]['EW'][iy, ix] = EW__ki[:, ii].mean() indices[name]['sigma_EW'][iy, ix] = EW__ki[:, ii].std() indices['SN'][iy, ix] = med_flux indices['e_SN'][iy, ix] = std_flux _a = [indices[k]['EW'][iy, ix] for k in k_list] _b = [indices[k]['sigma_EW'][iy, ix] for k in k_list] _a.append(indices['SN'][iy, ix]) _b.append(indices['e_SN'][iy, ix]) output_cube__Iyx[:, iy, ix] = np.append(_a, _b) return indices, output_cube__Iyx # TODO: # move to pyFIT3D/common/auto_ssp_tools.py # TODO: # create function auto_ssp_elines_rnd_cube_main() (structures version) # create function auto_ssp_elines_rnd_cube() (_main() helper for files) # TODO: # write documentation for the three new functions. def auto_ssp_desegmentation_main( wavelength, cube_flux, cube_eflux, output_files, ssp_file, config_file, out_file, seg_map, rss_nl_pars, rss_coeffs, ssp_nl_fit_file=None, sigma_inst=None, mask_list=None, elines_mask_file=None, w_min=None, w_max=None, nl_w_min=None, nl_w_max=None, R_V=None, extlaw=None, plot=None, fit_sigma_rnd=True, sps_class=StPopSynt, kdtree_nodes=4): """ """ # XXX: # All n_output_* should be defined from their "TO BE DEFINED" # dictionary sctructures (like _SSP_OUTPUT_INDEX) inside # pyFIT3D/common/constants.py n_output_spectra_list = 6 n_output_coeffs = 7 n_output_results = len(_SSP_OUTPUT_INDEX.keys()) redshift__s, sigma__s, AV__s = rss_nl_pars coeff__cs, norm__cs, e_norm__cs = rss_coeffs elines_out, coeffs_out, coeffs_nhlib_out, summary_out = output_files nc, ns = coeff__cs.shape nw, ny, nx = cube_flux.shape # Populate neighbors KDTree centcoor_yx_voxel__s = np.zeros((ns, 2), dtype='float') voxel_indices_yx__s = [] # kdtree_nodes means how many nodes the KDTree query will return, including the starting node (1). # I. e., each voxels is included in nearest_voxels__s list. k_nh = kdtree_nodes for id_voxel in range(1, ns + 1): yy, xx = np.where(seg_map == id_voxel) voxel_indices_yx__s.append([yy, xx]) centcoor_yx_voxel__s[id_voxel - 1] = [yy.mean(), xx.mean()] _, nearest_voxels__s = KDTree(centcoor_yx_voxel__s).query(centcoor_yx_voxel__s, k=kdtree_nodes) model_spectra = np.zeros((n_output_spectra_list, nw, ny, nx), dtype='float') results = np.zeros((n_output_results, ny, nx), dtype='float') results_coeffs = np.zeros((n_output_coeffs, nc, ny, nx), dtype='float') results_coeffs_nhlib = np.zeros((n_output_coeffs, kdtree_nodes, ny, nx), dtype='float') spx_nh_coeffs = np.zeros((ny, nx, kdtree_nodes, nc)) fitting_sequence = [] output_el_models = {} _tmpcf = ConfigAutoSSP(config_file) for i_s in range(_tmpcf.n_systems): system = _tmpcf.systems[i_s] elcf = _tmpcf.systems_config[i_s] k = f'{system["start_w"]}_{system["end_w"]}' output_el_models[k] = create_emission_lines_parameters(elcf, (ny, nx)) del _tmpcf # loop through voxels (CS bins) and its neiborhood (closest bins) # __i=0 for id_voxel in range(1, ns + 1): # if __i == 20: continue # __i += 1 i_v = id_voxel - 1 yy, xx = voxel_indices_yx__s[i_v] near_voxels = nearest_voxels__s[i_v] n_bins_voxel = near_voxels.shape ########################################################### # GUESSES AND COEFFS _redshift = [] _sigma = [] _AV = [] _coeffs = [] _sigma_coeffs = [] _min_coeffs = [] for _i_v in near_voxels: ## RETRIEVE GUESSES FROM CS AUTO SSP OUTPUT _redshift.append(redshift__s[_i_v]) _sigma.append(sigma__s[_i_v]) _AV.append(AV__s[_i_v]) _coeffs.append(coeff__cs[:, _i_v]) _sigma_coeffs.append(e_norm__cs[:, _i_v]) _min_coeffs.append(norm__cs[:, _i_v]) _redshift = np.asarray(_redshift) _sigma = np.asarray(_sigma) _AV = np.asarray(_AV) _coeffs = np.asarray(_coeffs).T _sigma_coeffs = np.asarray(_sigma_coeffs).T _min_coeffs = np.asarray(_min_coeffs).T _r_max = np.max(_redshift + 100/__c__) _r_min = np.min(_redshift - 100/__c__) _r_guess = [_redshift.mean(), 0.5*(_r_max-_r_min), _r_min, _r_max] _red, _d_red, _min_red, _max_red = _r_guess _s_max = 1.2*np.max(_sigma) _s_min = 0.8*np.min(_sigma) _s_guess = [_sigma.mean(), 1, _s_min, _s_max] _A_guess = [_AV.mean(), 0.15, 0, 1.6] ########################################################### ########################################################### # create temporary SSP from models! SSP_library = SSPModels(ssp_file) create_SSP_from_spectra(SSP_library, coeffs=_coeffs, output_filename='temp_ssp.fits') ########################################################### ########################################################### # Analyze spaxels for i_xy in zip(xx, yy): i_x, i_y = i_xy fitting_sequence.append(i_xy) f__w = copy(cube_flux[:, i_y, i_x]) ef__w = copy(cube_eflux[:, i_y, i_x]) spx_nh_coeffs[i_y, i_x] = _coeffs.T print(f"\n# ID {i_xy}/{(nx, ny)} ===============================================\n") cf, SPS = auto_ssp_elines_single_main( wavelength=wavelength, flux=f__w, eflux=ef__w, ssp_file='temp_ssp.fits', config_file=config_file, ssp_nl_fit_file=ssp_nl_fit_file, sigma_inst=sigma_inst, out_file=out_file, mask_list=mask_list, elines_mask_file=elines_mask_file, min=0-(0.33*f__w.mean()), max=1.05*f__w.max(), w_min=w_min, w_max=w_max, nl_w_min=nl_w_min, nl_w_max=nl_w_max, input_redshift=_red, delta_redshift=_d_red, min_redshift=_min_red, max_redshift=_max_red, input_sigma=_s_guess[0], delta_sigma=_s_guess[1], min_sigma=_s_guess[2], max_sigma=_s_guess[3], input_AV=_A_guess[0], delta_AV=_A_guess[1], min_AV=_A_guess[2], max_AV=_A_guess[3], plot=plot, single_ssp=False, spec_id=i_xy, fit_sigma_rnd=fit_sigma_rnd, sps_class=sps_class ) SPS.output_gas_emission(filename=elines_out, spec_id=i_xy) for system in SPS.config.systems: if system['EL'] is not None: k = f'{system["start_w"]}_{system["end_w"]}' append_emission_lines_parameters(system['EL'], output_el_models[k], i_xy) SPS.output_coeffs_MC(filename=coeffs_nhlib_out, write_header=id_voxel==1) # SPS.output(filename=summary_out, write_header=id_voxel==1, block_plot=False) model_spectra[:, :, i_y, i_x] = np.asarray(SPS.output_spectra_list) results_coeffs_nhlib[:, :, i_y, i_x] = np.asarray(SPS.output_coeffs) # results[:, i_y, i_x] = np.asarray(SPS.output_results) SPS.models = SSP_library SPS.ssp = SPS.models # _coeffs => _coeffs__cn (coeffs x neighb) _new_coeffs = np.dot(SPS.coeffs_norm, _coeffs.T) _new_min_coeffs = np.dot(SPS.min_coeffs_norm, _min_coeffs.T) _new_min_coeffs = _new_min_coeffs/_new_min_coeffs.sum() # new sigma coeffs 1 __x = SPS.coeffs_ssp_MC_rms[np.newaxis, :]*np.ones((SSP_library.n_models, k_nh)) _sc = np.sqrt((__x*_sigma_coeffs).sum(axis=1)) # new sigma coeffs 2 # _sc = np.asarray([np.sqrt(np.nansum((_c - np.nanmean(_c))**2)/(_c.size - 1)) for _c in _coeffs]) _new_sigma_coeffs = _sc SPS.coeffs_norm = _new_coeffs SPS.coeffs_ssp_MC_rms = _new_sigma_coeffs SPS.min_coeffs_norm = _new_min_coeffs SPS.coeffs_ssp_MC = _new_coeffs SPS.orig_best_coeffs = _new_min_coeffs # new_coeffs__cyx[:, i_y, i_x] = SPS.coeffs_norm # new_sigma_coeffs__cyx[:, i_y, i_x] = SPS.coeffs_ssp_MC_rms # _c_mass = SPS.coeffs_norm*SSP_library.mass_to_light # sum_mass = _c_mass.sum() # if sum_mass > 0: # cont_seg__yx[i_y, i_x] = 1 # new results SPS._MC_averages() wl_mass_factor = 3500 # TODO: # wl_mass_factor = w_max - w_min mass = SPS.mass_to_light*SPS.med_flux*wl_mass_factor lmass = 0 lml = 0 if mass: lmass = np.log10(mass) lml = np.log10(SPS.mass_to_light/wl_mass_factor) SPS.output_results = [ SPS.output_results[0], SPS.age_min, SPS.e_age_min, SPS.met_min, SPS.e_met_min, SPS.AV_min, SPS.e_AV_min, SPS.best_redshift, SPS.e_redshift, SPS.best_sigma, SPS.e_sigma, SPS.output_results[11], SPS.best_redshift, SPS.med_flux, SPS.rms, SPS.age_min_mass, SPS.e_age_min_mass, SPS.met_min_mass, SPS.e_met_min_mass, SPS.systemic_velocity, lml, lmass, ] SPS.output_coeffs = [ SPS.coeffs_norm, SPS.min_coeffs_norm, SPS.coeffs_ssp_MC, SPS.coeffs_ssp_MC_rms, SPS.coeffs_ssp_MC, SPS.coeffs_ssp_MC_rms, SPS.orig_best_coeffs, ] SPS.output_coeffs_MC(filename=coeffs_out, write_header=id_voxel==1) SPS.filename = ssp_file SPS.output(filename=summary_out, write_header=id_voxel==1, block_plot=False) results_coeffs[:, :, i_y, i_x] = np.asarray(SPS.output_coeffs) results[:, i_y, i_x] = np.asarray(SPS.output_results) return model_spectra, results, results_coeffs, results_coeffs_nhlib, spx_nh_coeffs, output_el_models, fitting_sequence def dump_auto_ssp_el_output(output_el_models, config_file, prefix): _tmpcf = ConfigAutoSSP(config_file) for i_s in range(_tmpcf.n_systems): system = _tmpcf.systems[i_s] elcf = _tmpcf.systems_config[i_s] k = f'{system["start_w"]}_{system["end_w"]}' _map_pref = f'map.auto_ssp.{prefix}.{k}.{name}' output_emission_lines_parameters(_map_pref, elcf, output_el_models[k]) del _tmpcf def map_auto_ssp_deseg(name, wavelength, model_spectra, results, results_coeffs, ssp_file, V__yx, mask_map__yx, prefix, output_prefix, wave_norm, inst_disp, sum_mass_age_out_file, cont_seg_file, scale_seg_file): _, nw, ny, nx = model_spectra.shape SSP_library = SSPModels(ssp_file) nc = SSP_library.n_models age__c = SSP_library.age_models met__c = SSP_library.metallicity_models # Starting clean maps! chi__yx = np.zeros((ny, nx)) age__yx = np.zeros((ny, nx)) age_mass__yx = np.zeros((ny, nx)) e_age__yx = np.zeros((ny, nx)) age_lum__yx = np.zeros((ny, nx)) e_age_lum__yx = np.zeros((ny, nx)) age_lum_M__yx = np.zeros((ny, nx)) e_age_lum_M__yx = np.zeros((ny, nx)) met_lum__yx = np.zeros((ny, nx)) met_mass__yx = np.zeros((ny, nx)) met_lumI__yx = np.zeros((ny, nx)) met_lumII__yx = np.zeros((ny, nx)) met_lumII_mass__yx = np.zeros((ny, nx)) e_met_lumII__yx = np.zeros((ny, nx)) met_e_lumII_mass__yx = np.zeros((ny, nx)) met__yx = np.zeros((ny, nx)) e_met__yx = np.zeros((ny, nx)) Av__yx = np.zeros((ny, nx)) e_Av__yx = np.zeros((ny, nx)) disp__yx = np.zeros((ny, nx)) vel__yx = np.zeros((ny, nx)) redshift__yx = np.zeros((ny, nx)) e_disp__yx = np.zeros((ny, nx)) e_vel__yx = np.zeros((ny, nx)) flux__yx = np.zeros((ny, nx)) e_flux__yx = np.zeros((ny, nx)) ml__yx = np.zeros((ny, nx)) mass__yx = np.zeros((ny, nx)) disp_km__yx = np.zeros((ny, nx)) e_disp_km__yx = np.zeros((ny, nx)) V_slice__yx = np.zeros((ny, nx)) scale__yx = np.zeros((ny, nx)) cont_seg__yx = np.zeros((ny, nx)) new_coeffs__cyx = np.zeros((nc, ny, nx)) new_sigma_coeffs__cyx = np.zeros((nc, ny, nx)) for ixy in itertools.product(range(nx),range(ny)): i_x, i_y = ixy raw_wave = wavelength model_ssp_min = model_spectra[1, :, i_y, i_x] output_results = results[:, i_y, i_x] output_coeffs = results_coeffs[:, :, i_y, i_x] # scale _sel_wave = trim_waves(raw_wave, [4500, 5500]) tmp_V_slice = model_ssp_min[_sel_wave] tmp_V_slice[~np.isfinite(tmp_V_slice)] = 0 V_slice__yx[i_y, i_x] = tmp_V_slice.mean() scale__yx[i_y, i_x] = V__yx[i_y, i_x]/V_slice__yx[i_y, i_x] _c_mass = output_results[0]*SSP_library.mass_to_light sum_mass = _c_mass.sum() if sum_mass > 0: cont_seg__yx[i_y, i_x] = 1 chi__yx[i_y, i_x] = output_results[_SSP_OUTPUT_INDEX['chi_joint']] age__yx[i_y, i_x] = output_results[_SSP_OUTPUT_INDEX['age_min']] met__yx[i_y, i_x] = output_results[_SSP_OUTPUT_INDEX['met_min']] Av__yx[i_y, i_x] = output_results[_SSP_OUTPUT_INDEX['AV_min']] redshift__yx[i_y, i_x] = output_results[_SSP_OUTPUT_INDEX['redshift']] _z = redshift__yx[i_y, i_x] _s = output_results[_SSP_OUTPUT_INDEX['sigma']] vel__yx[i_y, i_x] = _z*__c__ disp__yx[i_y, i_x] = _s/(1 + _z) if disp__yx[i_y, i_x] > inst_disp: disp_km__yx[i_y, i_x] = (np.sqrt(np.abs(disp__yx[i_y, i_x]**2 - inst_disp**2))/wave_norm)*__c__ else: disp_km__yx[i_y, i_x] = 0 flux__yx[i_y, i_x] = output_results[_SSP_OUTPUT_INDEX['med_flux']] e_age__yx[i_y, i_x] = output_results[_SSP_OUTPUT_INDEX['e_age_min']] e_met__yx[i_y, i_x] = output_results[_SSP_OUTPUT_INDEX['e_met_min']] e_Av__yx[i_y, i_x] = output_results[_SSP_OUTPUT_INDEX['e_AV_min']] e_vel__yx[i_y, i_x] = output_results[_SSP_OUTPUT_INDEX['e_redshift']]*__c__ e_disp__yx[i_y, i_x] = np.abs(output_results[_SSP_OUTPUT_INDEX['e_sigma']]) e_disp_km__yx[i_y, i_x] = (e_disp__yx[i_y, i_x]/wave_norm)*__c__ e_flux__yx[i_y, i_x] = output_results[_SSP_OUTPUT_INDEX['rms']] _iter = enumerate(zip(SSP_library.age_models, SSP_library.metallicity_models, np.log10(SSP_library.mass_to_light))) for i_c, (age, met, ml) in _iter: age_lum__yx[i_y, i_x] += output_coeffs[0][i_c]*(np.log10(age)) met_lum__yx[i_y, i_x] += output_coeffs[0][i_c]*met met_lumI__yx[i_y, i_x] += output_coeffs[1][i_c]*np.log10(met/__solar_metallicity__) met_lumII__yx[i_y, i_x] += output_coeffs[0][i_c]*np.log10(met/__solar_metallicity__) if sum_mass > 0: age_mass__yx[i_y, i_x] += _c_mass[i_c]*(np.log10(age)/sum_mass) ml__yx[i_y, i_x] += output_coeffs[0][i_c]*ml mass__yx[i_y, i_x] += output_coeffs[0][i_c]*(10**ml)*flux__yx[i_y, i_x]*1e-16 met_mass__yx[i_y, i_x] += _c_mass[i_c]*np.log10(met/__solar_metallicity__)/sum_mass if met__yx[i_y, i_x] > 0: e_met_lumII__yx[i_y, i_x] = e_met__yx[i_y, i_x]/met__yx[i_y, i_x] if age_lum__yx[i_y, i_x] > 0: e_age_lum__yx[i_y, i_x] = e_age__yx[i_y, i_x]/age_lum__yx[i_y, i_x] # from log Gyr to log yr age_lum__yx[i_y, i_x] = 9 + age_lum__yx[i_y, i_x] age_lum_M__yx[i_y, i_x] = 9 + age_mass__yx[i_y, i_x] met_lum__yx[i_y, i_x] = np.log10(met_lum__yx[i_y, i_x]/__solar_metallicity__) # Creating maps! print('Writting SFH CS fill maps...') chi_name = f'{output_prefix}_chi_ssp.fits.gz' array_to_fits(chi_name, chi__yx, overwrite=True) age_name = f'{output_prefix}_age_ssp.fits.gz' array_to_fits(age_name, age__yx, overwrite=True) age_name = f'{output_prefix}_age_mass_ssp.fits.gz' array_to_fits(age_name, age_mass__yx, overwrite=True) e_age_name = f'{output_prefix}_e_age_ssp.fits.gz' array_to_fits(e_age_name, e_age__yx, overwrite=True) age_name = f'{output_prefix}_log_age_yr_ssp.fits.gz' array_to_fits(age_name, age_lum__yx, overwrite=True) age_name = f'{output_prefix}_e_log_age_yr_ssp.fits.gz' array_to_fits(age_name, e_age_lum__yx, overwrite=True) age_name = f'{output_prefix}_log_age_yr_mass_ssp.fits.gz' array_to_fits(age_name, age_lum_M__yx, overwrite=True) met_name = f'{output_prefix}_met_ssp.fits.gz' array_to_fits(met_name, met__yx, overwrite=True) met_name = f'{output_prefix}_met_ZH_mass_ssp.fits.gz' array_to_fits(met_name, met_mass__yx, overwrite=True) met_name = f'{output_prefix}_met_ZH_ssp.fits.gz' array_to_fits(met_name, met_lumII__yx, overwrite=True) met_name = f'{output_prefix}_e_met_ZH_ssp.fits.gz' array_to_fits(met_name, e_met_lumII__yx, overwrite=True) Av_name = f'{output_prefix}_Av_ssp.fits.gz' array_to_fits(Av_name, Av__yx, overwrite=True) disp_name = f'{output_prefix}_disp_ssp.fits.gz' array_to_fits(disp_name, disp__yx, overwrite=True) disp_km_name = f'{output_prefix}_disp_km_h_ssp.fits.gz' array_to_fits(disp_km_name, disp_km__yx, overwrite=True) vel_name = f'{output_prefix}_vel_ssp.fits.gz' array_to_fits(vel_name, vel__yx, overwrite=True) flux_name = f'{output_prefix}_flux_ssp.fits.gz' array_to_fits(flux_name, flux__yx, overwrite=True) e_met_name = f'{output_prefix}_e_met_ssp.fits.gz' array_to_fits(e_met_name, e_met__yx, overwrite=True) e_Av_name = f'{output_prefix}_e_Av_ssp.fits.gz' array_to_fits(e_Av_name, e_Av__yx, overwrite=True) e_disp_name = f'{output_prefix}_e_disp_ssp.fits.gz' array_to_fits(e_disp_name, e_disp__yx, overwrite=True) e_disp_km_name = f'{output_prefix}_e_disp_km_h_ssp.fits.gz' array_to_fits(e_disp_km_name, e_disp_km__yx, overwrite=True) e_vel_name = f'{output_prefix}_e_vel_ssp.fits.gz' array_to_fits(e_vel_name, e_vel__yx, overwrite=True) e_flux_name = f'{output_prefix}_e_flux_ssp.fits.gz' array_to_fits(e_flux_name, e_flux__yx, overwrite=True) mass_name = f'{output_prefix}_mass_ssp.fits.gz' array_to_fits(mass_name, mass__yx, overwrite=True) ml_name = f'{output_prefix}_ml_ssp.fits.gz' array_to_fits(ml_name, ml__yx, overwrite=True) for i_c in range(nc): filename = f'{output_prefix}_NORM_{i_c}_ssp.fits.gz' sec__yx = new_coeffs__cyx[i_c] array_to_fits(filename, sec__yx, overwrite=True) write_img_header(filename, 'AGE', age__c[i_c]) write_img_header(filename, 'MET', met__c[i_c]) filename = f'{output_prefix}_eNORM_{i_c}_ssp.fits.gz' sec__yx = new_sigma_coeffs__cyx[i_c] array_to_fits(filename, sec__yx, overwrite=True) write_img_header(filename, 'AGE', age__c[i_c]) write_img_header(filename, 'MET', met__c[i_c]) output_prefix_crude = output_prefix.replace(name, 'NAME') age_unique = np.unique(SSP_library.age_models) n_age = age_unique.size met_unique = np.unique(SSP_library.metallicity_models) n_met = met_unique.size n_models = SSP_library.n_models _pref = output_prefix.replace('map.', '') f = open(f'sum_mass_age.{prefix}.{name}.pack.csv', 'w') f_out = open(sum_mass_age_out_file, 'w') _red = np.median(np.ma.masked_array(redshift__yx, mask=(redshift__yx != 0))) cosmo = get_cosmo(_red, __Hubble_constant__, __Omega_matter__, __Omega_Lambda__) ratio = 3.08567758e24 _luminosity_distance = 10**((cosmo['dm']-25)/5)*ratio _luminosity_factor = ((4*np.pi*(_luminosity_distance**2))*1e-16)/__solar_luminosity__ i_pack_out = 1 ML__yx = np.zeros_like(flux__yx) mass__yx = np.zeros_like(flux__yx) e_mass__yx = np.zeros_like(flux__yx) for im in range(n_models): age = SSP_library.age_models[im] met = SSP_library.metallicity_models[im] ML = SSP_library.mass_to_light[im] map__yx = get_data_from_fits(f'{output_prefix}_NORM_{im}_ssp.fits.gz') map__yx[~np.isfinite(map__yx)] = 0 e_map__yx = get_data_from_fits(f'{output_prefix}_eNORM_{im}_ssp.fits.gz') e_map__yx[~np.isfinite(e_map__yx)] = 0 _fact__yx = ML*flux__yx ML__yx += map__yx*ML mass__yx += map__yx*_fact__yx e_mass__yx += e_map__yx*_fact__yx if not im: cube__myx = np.zeros((n_models, ny, nx), dtype='float') e_cube__myx = np.zeros_like(cube__myx) cube__myx[im] = map__yx e_cube__myx[im] = e_map__yx print(f'{i_pack_out},{output_prefix_crude}_NORM_{im}_ssp.fits.gz,Luminosity Fraction for age-met {age:.4f}-{met:.4f} SSP, flux, fraction', file=f) i_pack_out += 1 cube__myx *= mask_map__yx e_cube__myx *= mask_map__yx Av__yx = get_data_from_fits(f'{output_prefix}_Av_ssp.fits.gz') lum__yx = flux__yx*_luminosity_factor*3500 mass__yx *= _luminosity_factor mass_to_light__yx = np.divide(mass__yx, lum__yx, where=lum__yx!=0, out=np.zeros((ny, nx))) e_mass__yx *= _luminosity_factor log10e = np.log10(np.exp(1)) loge10 = np.log(10) logMass__yx = np.log(mass__yx)/loge10 logMass__yx *= mask_map__yx array_to_fits(f'{output_prefix}_Mass_ssp.fits.gz', logMass__yx, overwrite=True) e_logMass__yx = np.divide(log10e*e_mass__yx, mass__yx, where=mass__yx!=0, out=np.zeros((ny, nx))) e_logMass__yx *= mask_map__yx array_to_fits(f'{output_prefix}_eMass_ssp.fits.gz', e_logMass__yx, overwrite=True) logMass_dust_corr__yx = logMass__yx + log10e*Av__yx array_to_fits(f'{output_prefix}_Mass_dust_cor_ssp.fits.gz', logMass_dust_corr__yx, overwrite=True) logML__yx = np.log(mass_to_light__yx)/loge10 array_to_fits(f'{output_prefix}_ML_ssp.fits.gz', logML__yx, overwrite=True) mass_sum = mass__yx.sum() e_mass_sum = e_mass__yx.sum() logMass_sum = np.log10(mass_sum) logeMass_sum = log10e*e_mass_sum/mass_sum logMass__yx[~np.isfinite(logMass__yx)] = 0 logMass_dust_corr__yx[~np.isfinite(logMass_dust_corr__yx)] = 0 mass_sum_dust_corr = np.where(np.abs(logMass_dust_corr__yx) > 0, 10**np.abs(logMass_dust_corr__yx), np.where(np.abs(logMass__yx) > 0, 10**np.abs(logMass__yx), 0)) logMass_dust_corr = np.log10(mass_sum_dust_corr.sum()) print(f'# name, log10(Mass)', file=f_out) print(f'Mass,{name},{logMass_sum},{logMass_dust_corr},{logeMass_sum}', file=f_out) for ia in range(n_age): age = age_unique[ia] age_t = f'{age:0>7.4f}' age__yx = np.zeros((ny, nx), dtype='float') e_age__yx = np.zeros((ny, nx), dtype='float') for im in range(n_models): if age == SSP_library.age_models[im]: age__yx += new_coeffs__cyx[im] e_age__yx += new_sigma_coeffs__cyx[im]**2 array_to_fits(f'{output_prefix}_{age_t}_NORM_age.fits.gz', age__yx, overwrite=True) array_to_fits(f'{output_prefix}_{age_t}_eNORM_age.fits.gz', np.sqrt(e_age__yx), overwrite=True) print(f'{i_pack_out},{output_prefix_crude}_{age_t}_NORM_age.fits.gz,Luminosity Fraction for age {age:.4f} SSP, flux, fraction', file=f) i_pack_out += 1 for iZ in range(n_met): met = met_unique[iZ] met_t = f'{met:0>6.4f}' met__yx = np.zeros((ny, nx), dtype='float') e_met__yx = np.zeros((ny, nx), dtype='float') for im in range(n_models): if met == SSP_library.metallicity_models[im]: met__yx += new_coeffs__cyx[im] e_met__yx += new_sigma_coeffs__cyx[im]**2 array_to_fits(f'{output_prefix}_{met_t}_NORM_met.fits.gz', met__yx, overwrite=True) array_to_fits(f'{output_prefix}_{met_t}_eNORM_met.fits.gz', np.sqrt(e_met__yx), overwrite=True) print(f'{i_pack_out},{output_prefix_crude}_{met_t}_NORM_met.fits.gz,Luminosity Fraction for met {met:.4f} SSP, flux, fraction', file=f) i_pack_out += 1 f_out.close() array_to_fits(f'{output_prefix}_V_ssp.fits.gz', V_slice__yx, overwrite=True) array_to_fits(cont_seg_file, cont_seg__yx, overwrite=True) array_to_fits(scale_seg_file, scale__yx, overwrite=True) def pack_dataproducts(header, files_dict, datacubes_names, output_filename): from pyFIT3D.common.constants import __version__ from astropy.io import fits import time from datetime import datetime SELECT_REG = fits.getdata(files_dict['SELECT_REG']) SELECT_REG_data = SELECT_REG/SELECT_REG hdul = fits.HDUList() # pyPipe3D signature header.append(card=('PIPELINE', __version__), end=True) # timestamp of datacube creation now = time.time() header.append(card=('UNIXTIME', int(now), f'{datetime.fromtimestamp(now)}'), end=True) hdul.append(fits.ImageHDU(data=None, header=header, name=datacubes_names[0])) for dcname in datacubes_names[1:]: fname = files_dict[dcname] if (fname is not None) & isfile(fname): t = fits.open(fname) _d = copy(t[0].data) _h = copy(t[0].header) if 'SSP' in dcname: try: i_ML = [_h[f'DESC_{i}'] for i in range(int(_h['NAXIS3']))].index(' average mass-to-light ratio of the stellar population') except ValueError: i_ML = 17 ML_med = np.nanmedian(_d[i_ML]) if ML_med > 2: _d[i_ML] -= np.log10(3500) if 'INDICES' in dcname: nI, ny, nx = _d.shape for _i in range(nI): _img = np.ma.masked_invalid(_d[_i], copy=True) _img *= SELECT_REG_data _img *= 1 _msk = _img == 0 _img = np.ma.masked_array(_img, mask=_msk) a = _img x, y = np.mgrid[0:a.shape[0], 0:a.shape[1]] _msk = np.ma.getmask(a) xygood = np.array((x[~_msk], y[~_msk])).T xybad = np.array((x[_msk], y[_msk])).T a[_msk] = a[~_msk][KDTree(xygood).query(xybad)[1]] _img = a _img *= SELECT_REG_data _img *= 1.0 _msk = _img == 0 _img = np.ma.masked_array(_img, mask=_msk, copy=True) _d[_i] = _img if ('FLUX_ELINES' in dcname) | ('SSP' in dcname) | ('SFH' in dcname) | ('INDICES' in dcname): _d[abs(_d) > 1e300] = np.nan hdul.append(fits.ImageHDU(data=_d, header=_h, name=dcname)) del(t) else: hdul.append(fits.ImageHDU(data=None, header=None, name=dcname)) hdul.writeto(output_filename, overwrite=True) def ana_single(name, input_config, force_z=None, x0=None, y0=None): #, only_binned_spaxels=False): config = input_config['general'] config_single = input_config['general_single_analysis'] config_seg = input_config['cube_segmentation'] config_elines = input_config['emission_lines'] config_pack = input_config['pack'] plot = config.getint('plot', 0) save_aux_files = config.getboolean('save_aux_files', True) dir_conf = config.get('dir_conf', '../legacy') dir_out = config.get('dir_out', '../out') dir_datacubes_out = config.get('dir_datacubes_out', '../out') inst_disp = config.getfloat('instrumental_dispersion') ########################################################################## # Load datacube # time_ini_loc = print_block_init('Loading datacube and calculating signal-to-noise...') cube_filename = config['cube_filename'].replace('NAME', name) org_cube__wyx, org_h, n_ext = get_data_from_fits(cube_filename, header=True, return_n_extensions=True) h_set_rss = {'CRVAL1': org_h['CRVAL3'], 'CDELT1': org_h['CDELT3'], 'CRPIX1': org_h['CRPIX3']} h_set_cube = {'CRVAL3': org_h['CRVAL3'], 'CDELT3': org_h['CDELT3'], 'CRPIX3': org_h['CRPIX3']} med_vel = org_h.get('MED_VEL') org_badpix__wyx = None org_error__wyx = None org_badspaxels__yx = None if org_h['EXTEND']: org_error__wyx = get_data_from_fits(cube_filename, extension=1) badpix_ext = config.getint('badpix_extension', None) if badpix_ext is not None: badpix_val = config.getint('badpix_value', 1) org_badpix__wyx = get_data_from_fits(cube_filename, extension=badpix_ext) if badpix_val != 1: _tmp_org_badpix__wyx = np.ones_like(org_badpix__wyx) _tmp_org_badpix__wyx[org_badpix__wyx != badpix_val] = 0 org_badpix__wyx = _tmp_org_badpix__wyx badpix_coll_threshold_frac = config.getfloat('badpix_coll_threshold_fraction', 0.5) badpix_outmask_value = config.getint('badpix_outmask_value', 1) org_badspaxels__yx = collapse_badpixels_mask(org_badpix__wyx, threshold_fraction=badpix_coll_threshold_frac, mask_value=badpix_outmask_value) else: print('# Unable to find a bad pixels extension on cube: ', flush=True) print('# - missing bad pixels spectra', flush=True) print('# - bad spaxels mask disable', flush=True) else: print('# Unable to find an EXTEND key on cube: ', flush=True) print('# - missing error spectra', flush=True) print('# - missing bad pixels spectra', flush=True) print('# - bad spaxels mask disable', flush=True) # original wavelength org_wave__w = get_wave_from_header(org_h, wave_axis=3) # Signal-to-noise information print('# -----=> get_SN_cube()', flush=True) SN__yx, signal__yx, noise__yx = get_SN_cube(copy(org_cube__wyx), org_wave__w) if save_aux_files: with open(f'SYS_VEL_NOW_{name}', 'w') as f: print(f'{cube_filename} {med_vel}', file=f) SN_map_filename = f'map_SN.{name}.fits.gz' signal_map_filename = f'signal.{name}.fits.gz' noise_map_filename = f'noise.{name}.fits.gz' array_to_fits(SN_map_filename, SN__yx, overwrite=True) array_to_fits(signal_map_filename, signal__yx, overwrite=True) array_to_fits(noise_map_filename, noise__yx, overwrite=True) print_done(time_ini=time_ini_loc) ########################################################################## ########################################################################## # pseudo V-band slice time_ini_loc = print_block_init('Creating pseudo V-band map...') slice_prefix = f'img_{name}' slice_config = config.get('slice_config_file') # get_slice(cube_filename, slice_prefix, slice_config) print('# -----=> get_slice()', flush=True) slices = get_slice(copy(org_cube__wyx), org_wave__w, slice_prefix, slice_config) V_img_slice_key = f'{slice_prefix}_V_4500_5500' #.fits' if save_aux_files: V_slice_filename = f'{slice_prefix}_V_4500_5500.fits' array_to_fits(V_slice_filename, slices[V_img_slice_key], overwrite=True) print_done(time_ini=time_ini_loc) ########################################################################## ########################################################################## # Creating mask # V__yx = copy(slices[V_img_slice_key]) ny, nx = V__yx.shape nx1, ny1 = nx, ny # V__yx = get_data_from_fits(tmp_img_filename) msk__yx = np.zeros_like(V__yx) # EL: deprecated - MaNGA v1.4 mask msk_filename = config.get('spatial_mask_file', None) if msk_filename is not None: msk_filename = msk_filename.replace('NAME', name) if isfile(msk_filename): msk__yx = get_data_from_fits(msk_filename) ny1, nx1 = msk__yx.shape print(f'# Dimensions = {nx},{ny},{nx1},{ny1}') # mask_val = 1 means that this is not used mask_stars = config.getboolean('mask_stars', False) #config.getint('mask_stars', 1) mask_val = 1 if mask_stars: mask_val = 0 V_msk__yx = mask_stars_in_FoV(name, copy(V__yx), msk__yx, mask_val=mask_val) if save_aux_files: msk_filename = f'{name}.mask.fits' array_to_fits(msk_filename, V_msk__yx, overwrite=True) # Appending badpix mask V_msk__yx: if some spaxel has only bad pixels, # this spaxel should be masked. if org_badspaxels__yx is not None: V_msk__yx[org_badspaxels__yx == 1] = 0 if save_aux_files: msk_filename = f'{name}.mask.nobadspaxels.fits' array_to_fits(msk_filename, V_msk__yx, overwrite=True) ########################################################################## ########################################################################## # Central spectrum integrated spectrum # # pseudo V-band map time_ini_loc = print_block_init('Analyzing central and integrated spectra...') if x0 is None: x0 = config_single.getfloat('x0', None) if y0 is None: y0 = config_single.getfloat('y0', None) if x0 is None or y0 is None: nx0 = int(nx*.33) nx1 = int(nx*2*.33) ny0 = int(ny*.33) ny1 = int(ny*2*.33) iy, ix = np.indices(V__yx.shape) iy = iy[ny0:ny1, nx0:nx1] ix = ix[ny0:ny1, nx0:nx1] sel = V_msk__yx[ny0:ny1, nx0:nx1] != 2 V_slice = copy(V__yx[ny0:ny1, nx0:nx1][sel]) __V = V_slice**5 _s = __V.sum() x0 = ix[sel].dot(__V)/_s y0 = iy[sel].dot(__V)/_s # clean_nan(V_img_filename, badval=-1) V__yx[np.isnan(V__yx)] = -1 print(f'# x0:{x0} y0:{y0}') # creating spectra (central and integrated) print('# -----=> radial_sum_cube_e()', flush=True) # radial_sum_cube_e(cube_filename, delta_R=2.5, x0=x0, y0=y0, output_rss_fits=radial_sum_cube_filename) # (EADL) To think: # radial_sum_cube_e() do not consider V-band mask (which inlcudes masked # stars in the field), only a bad pixels mask (input_mask). The spatial mask # should also be processed by radial_sum_cube_e() function. r = radial_sum_cube_e(copy(org_cube__wyx), x0=x0, y0=y0, delta_R=config_single.getfloat('radial_sum_cube_delta_R', 2.5), input_mask=org_badpix__wyx, input_error=org_error__wyx) output__rw, e_output__rw, mask_output__rw = r print('# -----=> img2spec_e() central spectrum', flush=True) # img2spec_e(radial_sum_cube_filename, 0, cen_spec_filename) f_cen__w = copy(output__rw[0]) f_cen__w[~np.isfinite(f_cen__w)] = 0 sqrt_f_cen__w = np.sqrt(np.abs(f_cen__w)) e_f_cen__w = copy(e_output__rw[0]) e_f_cen__w[~np.isfinite(e_f_cen__w)] = 1e12 e_f_cen__w = np.where(e_f_cen__w > sqrt_f_cen__w, sqrt_f_cen__w, e_f_cen__w) print('# -----=> img2spec_e() integrated spectrum', flush=True) # img2spec_e(radial_sum_cube_filename, 5, int_spec_filename) f_int__w = copy(output__rw[5]) f_int__w[~np.isfinite(f_int__w)] = 0 sqrt_f_int__w = np.sqrt(np.abs(f_int__w)) e_f_int__w = copy(e_output__rw[5]) e_f_int__w[~np.isfinite(e_f_int__w)] = 1e12 e_f_int__w = np.where(e_f_int__w > sqrt_f_int__w, sqrt_f_int__w, e_f_int__w) # Needed final files cen_spec_filename = config_single.get('cen_spec_filename', None) cen_spec_filename = f'{name}.spec_5.txt' if cen_spec_filename is None else cen_spec_filename.replace('NAME', name) output_spectra(org_wave__w, [f_cen__w, e_f_cen__w], cen_spec_filename) int_spec_filename = config_single.get('int_spec_filename', None) int_spec_filename = f'{name}.spec_30.txt' if int_spec_filename is None else int_spec_filename.replace('NAME', name) output_spectra(org_wave__w, [f_int__w, e_f_int__w], int_spec_filename) if save_aux_files: radial_sum_cube_filename = f'rad.{name}.rss.fits' e_radial_sum_cube_filename = f'e_{radial_sum_cube_filename}' weight_radial_sum_cube_filename = f'weight.{radial_sum_cube_filename}' array_to_fits(radial_sum_cube_filename, output__rw, header=h_set_rss, overwrite=True) array_to_fits(e_radial_sum_cube_filename, e_output__rw, header=h_set_rss, overwrite=True) array_to_fits(weight_radial_sum_cube_filename, mask_output__rw, header=h_set_rss, overwrite=True) # CLEAN MEMORY del(r) del(output__rw) del(e_output__rw) del(mask_output__rw) gc.collect() # define redshift range ############################################################################ # (EADL): To think: This entire block has not been used. ############################################################################ ############################################################################ # print('# -----=> vel_eline_spec()', flush=True) # # npeaks, vel = vel_eline_spec(cen_spec_filename, 2, 0.7, 'junk', 3728, 3800, 4000, output_file=f'vel_eline_spec.{name}.txt') # w_min = config_single.getint('vel_eline_min_wavelength', 3800) # w_max = config_single.getint('vel_eline_max_wavelength', 4000) # wave_ref = config_single.getfloat('vel_eline_reference_wavelength', 3728) # nsearch = config_single.getint('vel_eline_nsearch', 2) # imin = config_single.getfloat('vel_eline_imin', 0.7) # sel__w = trim_waves(org_wave__w, [w_min, w_max]) # vel, _, npeaks = vel_eline(f_cen__w[sel__w], org_wave__w[sel__w], nsearch=nsearch, imin=imin, wave_ref=wave_ref) # if npeaks == 1: # red = vel/__c__ # red_f = 300/__c__ # min_red = red - red_f # min_red = red + red_f # else: # # EL: not used in MaNGA I supose, because else is just rewritten # # below, in NSA redshift if. # red = 0.02 # min_red = 0.0045 # max_red = 0.032 ############################################################################ if force_z is None: force_z = config_single.getfloat('force_redshift', None) # (EADL) TODO: Create a parser to read a redshift file (^NAME\tXXXX$ probable format) # if force_z is None: # force_z = get_redshift_from_file(config_single.get('redshift_file', None)) # print(force_z) if (force_z is not None) and (force_z > 0): red = force_z red_hr_kms = config_single.getfloat('redshift_half_range', 600) red_f = red_hr_kms/__c__ min_red = red - red_f max_red = red + red_f d_red = 0.2*red_hr_kms/__c__ print(f'# USING FORCED REDSHIFT') else: red = config_single.getfloat('input_redshift', 0.02) min_red = config_single.getfloat('min_redshift', -0.001) max_red = config_single.getfloat('max_redshift', 0.2) d_red = config_single.getfloat('delta_redshift_kms', 300)/__c__ print(f'# redshift = {red:.8f} - min = {min_red:.8f} - max = {max_red:.8f} - delta = {d_red:.8f}') time_ini_auto_ssp = print_block_init('auto_ssp_elines_rnd(): Analyzing central spectrum...') tmp_cf = input_config['cen_auto_ssp_no_lines'] auto_ssp_output_filename = tmp_cf.get('output_filename', 'auto_ssp_Z.NAME.cen.out').replace('NAME', name) remove_isfile(auto_ssp_output_filename) auto_ssp_output_elines_filename = 'elines_' + auto_ssp_output_filename remove_isfile(auto_ssp_output_elines_filename) auto_ssp_output_coeffs_filename = 'coeffs_' + auto_ssp_output_filename remove_isfile(auto_ssp_output_coeffs_filename) auto_ssp_output_fits_filename = 'output.' + auto_ssp_output_filename + '.fits' remove_isfile(auto_ssp_output_fits_filename) w_min = tmp_cf.getint('min_wavelength', 3800) w_max = tmp_cf.getint('max_wavelength', 7000) nl_w_min = tmp_cf.getint('nonlinear_min_wavelength', 3850) nl_w_max = tmp_cf.getint('nonlinear_max_wavelength', 4700) cf, SPS = auto_ssp_elines_single_main( copy(org_wave__w), copy(f_cen__w), copy(e_f_cen__w), tmp_cf.get('models_file'), tmp_cf.get('config_file'), auto_ssp_output_filename, ssp_nl_fit_file=tmp_cf.get('nonlinear_models_file'), mask_list=tmp_cf.get('mask_file', None), elines_mask_file=tmp_cf.get('emission_lines_mask_file', None), min=None, max=None, w_min=w_min, w_max=w_max, nl_w_min=nl_w_min, nl_w_max=nl_w_max, input_redshift=red, delta_redshift=d_red, min_redshift=min_red, max_redshift=max_red, input_sigma=tmp_cf.getfloat('input_sigma'), delta_sigma=tmp_cf.getfloat('delta_sigma'), min_sigma=tmp_cf.getfloat('min_sigma'), max_sigma=tmp_cf.getfloat('max_sigma'), input_AV=tmp_cf.getfloat('input_AV'), delta_AV=tmp_cf.getfloat('delta_AV'), min_AV=tmp_cf.getfloat('min_AV'), max_AV=tmp_cf.getfloat('max_AV'), R_V=tmp_cf.getfloat('RV'), extlaw=tmp_cf.get('extlaw'), plot=plot, fit_sigma_rnd=tmp_cf.getboolean('fit_sigma_rnd', True), fit_gas=False, ) SPS.output_gas_emission(filename=auto_ssp_output_elines_filename) SPS.output_fits(filename=auto_ssp_output_fits_filename) SPS.output_coeffs_MC(filename=auto_ssp_output_coeffs_filename) SPS.output(filename=auto_ssp_output_filename, block_plot=False) print_done(time_ini=time_ini_auto_ssp, message='auto_ssp_elines_rnd(): DONE!') # data = np.genfromtxt(auto_ssp_output_filename, usecols=[7,9], delimiter=',', # dtype=[('redshift', 'float'), ('disp', 'float')]) # red = data['redshift'] # disp_max = 1.5*data['disp'] # read redshift and dispersion from auto_ssp output red = SPS.best_redshift disp_max = 1.5*SPS.best_sigma # ################################################################### # # A faster way to derive the redshift ############################# # ################################################################### # # This entire block can replace the previous # # auto_ssp_elines_single_main() call. # ################################## # # Do not perform SPS, only fits redshift and sigma # # also do not produce output files # ################################################################### # from pyFIT3D.modelling.stellar import StPopSynt # tmp_cf = input_config['cen_auto_ssp_no_lines'] # sigma_set = [ # tmp_cf.getfloat('input_sigma'), tmp_cf.getfloat('delta_sigma'), # tmp_cf.getfloat('min_sigma'), tmp_cf.getfloat('max_sigma') # ] # redshift_set = [red, d_red, min_red, max_red] # # do not fit AV # AV_set = [tmp_cf.getfloat('input_AV'), 0, 0, 0] # cf = ConfigAutoSSP(tmp_cf.get('config_file'), redshift_set=redshift_set, sigma_set=sigma_set, AV_set=AV_set) # w_min = tmp_cf.getint('min_wavelength', 3800) # w_max = tmp_cf.getint('max_wavelength', 7000) # nl_w_min = tmp_cf.getint('nonlinear_min_wavelength', 3850) # nl_w_max = tmp_cf.getint('nonlinear_max_wavelength', 4700) # ssp_nl_fit_file=, # mask_list=tmp_cf.get('mask_file', None), # elines_mask_file=tmp_cf.get('emission_lines_mask_file', None), # SPS = StPopSynt(config=cf, wavelength=copy(org_wave__w), # flux=copy(f_cen__w), eflux=copy(e_f_cen__w), # mask_list=tmp_cf.get('mask_file', None), # elines_mask_file=tmp_cf.get('emission_lines_mask_file', None), # ssp_file=tmp_cf.get('models_file'), # ssp_nl_fit_file=tmp_cf.get('nonlinear_models_file'), # w_min=w_min, w_max=w_max, nl_w_min=nl_w_min, nl_w_max=nl_w_max, # R_V=tmp_cf.getfloat('RV'), extlaw=tmp_cf.get('extlaw'), # sigma_inst=None, out_file=None, fit_gas=False, plot=plot) # msg_cut = f' - cut value: {cf.CUT_MEDIAN_FLUX:6.4f}' # if cf.CUT_MEDIAN_FLUX == 0: # msg_cut = ' - Warning: no cut (CUT_MEDIAN_FLUX = 0)' # print(f'-> median raw flux = {SPS.median_flux:6.4f}{msg_cut}') # SPS.cut = False # if SPS.median_flux > cf.CUT_MEDIAN_FLUX: # and (median_flux > cf.ABS_MIN): # # redshift # SPS.non_linear_fit(is_guided_sigma, fit_sigma_rnd=fit_sigma_rnd) # red = SPS.best_redshift # disp_max = 1.5*SPS.best_sigma # else: # # max sigma is in AA -> km/s using 5000 A as reference wl # disp_max = tmp_cf.getfloat('max_sigma')*5000/__c__ # ################################################################### red_f = 300/__c__ min_red = red - red_f max_red = red + red_f d_red = 0.1*red_f vel = red*__c__ print(f'# Redshift = {red:.8f} - Vel = {vel:.8f}') # fix redshift at config file config_file = f'auto.{name}.config' print('# -----=> redshift_config()', flush=True) redshift_config( input_config=input_config['auto_ssp_config'].get('strong'), redshift=red, output_config=config_file ) # Analyze central spectrum with no mask time_ini_auto_ssp = print_block_init('auto_ssp_elines_rnd(): Analyzing central spectrum with instrumental sigma (no mask)...') tmp_cf = input_config['cen_auto_ssp_inst_disp'] w_min = w_min*(1 + red) w_max = w_max*(1 + red) nl_w_min = nl_w_min*(1 + red) nl_w_max = nl_w_max*(1 + red) org_w_minmax = [w_min, w_max] org_nl_w_minmax = [w_min, w_max] auto_ssp_output_filename = tmp_cf.get('no_mask_output_filename', 'auto_ssp_no_mask.NAME.cen.out').replace('NAME', name) remove_isfile(auto_ssp_output_filename) auto_ssp_output_elines_filename = 'elines_' + auto_ssp_output_filename remove_isfile(auto_ssp_output_elines_filename) auto_ssp_output_coeffs_filename = 'coeffs_' + auto_ssp_output_filename remove_isfile(auto_ssp_output_coeffs_filename) auto_ssp_output_fits_filename = 'output.' + auto_ssp_output_filename + '.fits' remove_isfile(auto_ssp_output_fits_filename) cen_auto_ssp_output_fits_filename = auto_ssp_output_fits_filename cf, SPS = auto_ssp_elines_single_main( wavelength=copy(org_wave__w), flux=copy(f_cen__w), eflux=copy(e_f_cen__w), ssp_file=tmp_cf.get('models_file'), ssp_nl_fit_file=tmp_cf.get('nonlinear_models_file'), config_file=config_file, out_file=auto_ssp_output_filename, sigma_inst=tmp_cf.getfloat('instrumental_dispersion'), mask_list=None, elines_mask_file=tmp_cf.get('emission_lines_mask_file', None), min=-3, max=50, w_min=w_min, w_max=w_max, nl_w_min=nl_w_min, nl_w_max=nl_w_max, input_redshift=red, delta_redshift=0, min_redshift=min_red, max_redshift=max_red, input_sigma=tmp_cf.getfloat('input_sigma'), delta_sigma=tmp_cf.getfloat('delta_sigma'), min_sigma=tmp_cf.getfloat('min_sigma'), max_sigma=tmp_cf.getfloat('max_sigma'), input_AV=tmp_cf.getfloat('input_AV'), delta_AV=tmp_cf.getfloat('delta_AV'), min_AV=tmp_cf.getfloat('min_AV'), max_AV=tmp_cf.getfloat('max_AV'), R_V=tmp_cf.getfloat('RV'), extlaw=tmp_cf.get('extlaw'), plot=plot, fit_sigma_rnd=tmp_cf.getboolean('fit_sigma_rnd', True), ) SPS.output_gas_emission(filename=auto_ssp_output_elines_filename) SPS.output_fits(filename=auto_ssp_output_fits_filename) SPS.output_coeffs_MC(filename=auto_ssp_output_coeffs_filename) SPS.output(filename=auto_ssp_output_filename, block_plot=False) print_done(time_ini=time_ini_auto_ssp, message='auto_ssp_elines_rnd(): DONE!') # Analyze central spectrum with mask time_ini_auto_ssp = print_block_init('auto_ssp_elines_rnd(): Analyzing central spectrum with instrumental sigma...') auto_ssp_output_filename = tmp_cf.get('output_filename', 'auto_ssp.NAME.cen.out').replace('NAME', name) remove_isfile(auto_ssp_output_filename) auto_ssp_output_elines_filename = 'elines_' + auto_ssp_output_filename remove_isfile(auto_ssp_output_elines_filename) auto_ssp_output_coeffs_filename = 'coeffs_' + auto_ssp_output_filename remove_isfile(auto_ssp_output_coeffs_filename) auto_ssp_output_fits_filename = 'output.' + auto_ssp_output_filename + '.fits' remove_isfile(auto_ssp_output_fits_filename) cf, SPS = auto_ssp_elines_single_main( wavelength=copy(org_wave__w), flux=copy(f_cen__w), eflux=copy(e_f_cen__w), ssp_file=tmp_cf.get('models_file'), ssp_nl_fit_file=tmp_cf.get('nonlinear_models_file'), config_file=config_file, out_file=auto_ssp_output_filename, sigma_inst=tmp_cf.getfloat('instrumental_dispersion'), mask_list=tmp_cf.get('mask_file', None), elines_mask_file=tmp_cf.get('emission_lines_mask_file', None), min=-3, max=50, w_min=w_min, w_max=w_max, nl_w_min=nl_w_min, nl_w_max=nl_w_max, input_redshift=red, delta_redshift=0, min_redshift=min_red, max_redshift=max_red, input_sigma=tmp_cf.getfloat('input_sigma'), delta_sigma=tmp_cf.getfloat('delta_sigma'), min_sigma=tmp_cf.getfloat('min_sigma'), max_sigma=tmp_cf.getfloat('max_sigma'), input_AV=tmp_cf.getfloat('input_AV'), delta_AV=tmp_cf.getfloat('delta_AV'), min_AV=tmp_cf.getfloat('min_AV'), max_AV=tmp_cf.getfloat('max_AV'), R_V=tmp_cf.getfloat('RV'), extlaw=tmp_cf.get('extlaw'), plot=plot, fit_sigma_rnd=tmp_cf.getboolean('fit_sigma_rnd', True), ) SPS.output_gas_emission(filename=auto_ssp_output_elines_filename) SPS.output_fits(filename=auto_ssp_output_fits_filename) SPS.output_coeffs_MC(filename=auto_ssp_output_coeffs_filename) SPS.output(filename=auto_ssp_output_filename, block_plot=False) print_done(time_ini=time_ini_auto_ssp, message='auto_ssp_elines_rnd(): DONE!') # Analyze integrated spectrum # data = np.genfromtxt(auto_ssp_output_filename, usecols=[9], delimiter=',', # dtype=[('disp', 'float')]) time_ini_auto_ssp = print_block_init('auto_ssp_elines_rnd(): Analyzing integrated spectrum with instrumental sigma...') tmp_cf = input_config['int_auto_ssp_inst_disp'] sigma = SPS.best_sigma max_sigma = sigma + 30 guess_sigma = sigma print(f'# Max sigma:{max_sigma:.8f} - Guess sigma:{guess_sigma:.8f}') guess_sigma = 30 auto_ssp_output_filename = tmp_cf.get('output_filename', 'auto_ssp.NAME.int.out').replace('NAME', name) remove_isfile(auto_ssp_output_filename) int_auto_ssp_output_filename = auto_ssp_output_filename auto_ssp_output_elines_filename = 'elines_' + auto_ssp_output_filename remove_isfile(auto_ssp_output_elines_filename) auto_ssp_output_coeffs_filename = 'coeffs_' + auto_ssp_output_filename remove_isfile(auto_ssp_output_coeffs_filename) auto_ssp_output_fits_filename = 'output.' + auto_ssp_output_filename + '.fits' remove_isfile(auto_ssp_output_fits_filename) cf, SPS = auto_ssp_elines_single_main( wavelength=copy(org_wave__w), flux=copy(f_int__w), eflux=copy(e_f_int__w), ssp_file=tmp_cf.get('models_file'), ssp_nl_fit_file=tmp_cf.get('nonlinear_models_file'), config_file=config_file, out_file=auto_ssp_output_filename, sigma_inst=tmp_cf.getfloat('instrumental_dispersion'), mask_list=tmp_cf.get('mask_file', None), elines_mask_file=tmp_cf.get('emission_lines_mask_file', None), min=-3, max=50, w_min=w_min, w_max=w_max, nl_w_min=nl_w_min, nl_w_max=nl_w_max, input_redshift=red, delta_redshift=d_red, min_redshift=min_red, max_redshift=max_red, input_sigma=guess_sigma, delta_sigma=tmp_cf.getfloat('delta_sigma'), min_sigma=tmp_cf.getfloat('min_sigma'), max_sigma=max_sigma, input_AV=tmp_cf.getfloat('input_AV'), delta_AV=tmp_cf.getfloat('delta_AV'), min_AV=tmp_cf.getfloat('min_AV'), max_AV=tmp_cf.getfloat('max_AV'), R_V=tmp_cf.getfloat('RV'), extlaw=tmp_cf.get('extlaw'), plot=plot, fit_sigma_rnd=tmp_cf.getboolean('fit_sigma_rnd', True), ) SPS.output_gas_emission(filename=auto_ssp_output_elines_filename) SPS.output_fits(filename=auto_ssp_output_fits_filename) SPS.output_coeffs_MC(filename=auto_ssp_output_coeffs_filename) SPS.output(filename=auto_ssp_output_filename, block_plot=False) print_done(time_ini=time_ini_auto_ssp, message='auto_ssp_elines_rnd(): DONE!') print_done(time_ini=time_ini_loc, message='Central and integrated analysis DONE!') ########################################################################## ########################################################################## # FoV segmentation (RSS build) # time_ini_loc = print_block_init('FoV segmentation (RSS build)...') print('# -----=> med2df()', flush=True) # med2df(input_fits=V_img_filename, # x_width=2, y_width=2, # output_fits=filt_V_img_filename) # imarith(filt_V_img_filename, '*', msk_filename, filt_V_img_filename) filt_V__yx = median_filter(V__yx, size=(2, 2), mode='reflect') mV__yx = filt_V__yx*V_msk__yx ######## # (EADL) # NOTE: config[general_single_analysis][pseudo_V_band_file] usage is # deprecated. _tmp_V_img_filename = config_single.get('pseudo_V_band_file', 'NAME.V.fits') ######## V_img_filename = config_seg.get('pseudo_V_band_file', _tmp_V_img_filename).replace('NAME', name) array_to_fits(V_img_filename, mV__yx, overwrite=True) # Segmentation map from signal and noise maps print('# -----=> cont_seg_all_SN()', flush=True) flux_limit_peak = config_seg.getfloat('flux_limit_peak', 0.001) target_SN = config_seg.getfloat('target_SN', 50) min_SN = config_seg.getfloat('min_SN', 1) frac_peak = config_seg.getfloat('frac_peak', 0.85) min_flux = config_seg.getfloat('min_flux', 0.001) seg_map__yx, dmask_map__yx = cont_seg_all_SN(copy(signal__yx), copy(noise__yx), flux_limit_peak=flux_limit_peak, target_SN=target_SN, min_SN=min_SN, frac_peak=frac_peak, min_flux=min_flux) # Do not get stucked due to bad data. if (seg_map__yx.sum() == 0) or (dmask_map__yx.sum() == (nx*ny)): print(f'{basename(sys.argv[0])}: cont_seg_all_SN: no segmentation. Stopping analysis.') return # mask based in SN < 1 mask_map__yx = (SN__yx > min_SN).astype(int) # extract segmented spectra print('# -----=> spec_extract_cube_mean()', flush=True) # spec_extract_cube_mean(cube_filename, cont_seg_filename, cont_seg_rss_filename) org_rss__sw, _, _x_mean, _y_mean, npt_mean = spec_extract_cube_mean(copy(org_cube__wyx), seg_map__yx, badpix__wyx=org_badpix__wyx) ns = org_rss__sw.shape[0] print('# -----=> spec_extract_cube_error()', flush=True) # spec_extract_cube_error(cube_filename, cont_seg_filename, e_cont_seg_rss_filename) org_e_rss__sw, inv_seg__yx, _x_error, _y_error, npt_error = spec_extract_cube_error(copy(org_error__wyx), seg_map__yx, badpix__wyx=org_badpix__wyx, fov_selection__yx=mV__yx.astype('bool')) # org_rss__sw, h_org_rss = get_data_from_fits(cont_seg_rss_filename, header=True) # wave_rss__w = get_wave_from_header(h_org_rss, wave_axis=1) wave_rss__w = copy(org_wave__w) # get signal to noise from RSS print('# -----=> get_SN_rss()', flush=True) SN__s, signal__s, noise__s = get_SN_rss(copy(org_rss__sw), wave_rss__w) # output SN rss file # CSV to FITS print('# -----=> csv_to_map_seg()', flush=True) # csv_to_map_seg(output_SN_rss_csv_filename, 1, cont_seg_filename, output_SN_rss_fits_filename) SN_rss__yx, norm_SN_rss__yx, area_rss__yx = SN_map_seg(copy(SN__s), seg_map__yx) # RSS SEG to Cube print('# -----=> rss_seg2cube()', flush=True) cube__wyx = rss_seg2cube(org_rss__sw, seg_map__yx) # get pseudo V-band slice print('# -----=> get_slice()', flush=True) slice_prefix = f'SEG_img_{name}' slices = get_slice(copy(cube__wyx), wave_rss__w, slice_prefix, slice_config) k = list(slices.keys())[0] V_slice__yx = slices[k] V_slice__yx[np.isnan(V_slice__yx)] = -1 scale_seg__yx = mV__yx/V_slice__yx # Needed final files cont_seg_filename = config_seg.get('cont_seg_file', f'cont_seg.NAME.fits').replace('NAME', name) array_to_fits(cont_seg_filename, seg_map__yx, overwrite=True) scale_seg_filename = config_seg.get('scale_seg_file', f'scale.seg.NAME.fits').replace('NAME', name) array_to_fits(scale_seg_filename, scale_seg__yx, overwrite=True) SN_mask_map_filename = f'mask.{name}.V.fits' array_to_fits(SN_mask_map_filename, mask_map__yx, overwrite=True) if save_aux_files: filt_V_img_filename = f'm{V_img_filename}' array_to_fits(filt_V_img_filename, mV__yx, overwrite=True) diffuse_seg_filename = f'DMASK.{name}.fits' array_to_fits(diffuse_seg_filename, dmask_map__yx, overwrite=True) cont_seg_rss_filename = f'CS.{name}.RSS.fits' array_to_fits(cont_seg_rss_filename, org_rss__sw, header=h_set_rss, overwrite=True) size = np.sqrt(nx**2+ny**2)/(2*ns) output_rss_txt = cont_seg_rss_filename.replace('fits', 'pt.txt') with open(output_rss_txt, 'w') as f: output_header = f'C {size} {size} 0' print(output_header, file=f) # output_header = f'(1) id\n(2) S/N\n(3) Signal\n(4) Noise' np.savetxt(f, list(zip(list(range(ns)), _x_mean/npt_mean, _y_mean/npt_mean, [1]*ns)), fmt=['%d'] + 2*['%.18g'] + ['%d'], delimiter=' ') print(f'# {cont_seg_rss_filename} and {output_rss_txt} created') e_cont_seg_rss_filename = f'e_{cont_seg_rss_filename}' array_to_fits(e_cont_seg_rss_filename, org_e_rss__sw, header=h_set_rss, overwrite=True) size = np.sqrt(nx**2+ny**2)/(2*ns) output_rss_txt = e_cont_seg_rss_filename.replace('fits', 'pt.txt') with open(output_rss_txt, 'w') as f: output_header = f'C {size} {size} 0' print(output_header, file=f) # output_header = f'(1) id\n(2) S/N\n(3) Signal\n(4) Noise' np.savetxt(f, list(zip(list(range(ns)), _x_error/npt_error, _y_error/npt_error, [1]*ns)), fmt=['%d'] + 2*['%.18g'] + ['%d'], delimiter=' ') print(f'# {e_cont_seg_rss_filename} and {output_rss_txt} created') output_SN_rss_csv_filename = f'SN_{name}.CS.rss.csv' with open(output_SN_rss_csv_filename, 'w') as f: output_header = f'(1) id\n(2) S/N\n(3) Signal\n(4) Noise' np.savetxt(f, list(zip(list(range(SN__s.size)), SN__s, signal__s, noise__s)), fmt=['%d'] + 3*['%.18g'], header=output_header, delimiter=',') output_SN_rss_fits_filename = f'SN_{name}.CS.fits' output_norm_SN_rss_fits_filename = f'norm_{output_SN_rss_fits_filename}' output_area_rss_fits_filename = f'area_{output_SN_rss_fits_filename}' array_to_fits(output_SN_rss_fits_filename, SN_rss__yx, overwrite=True) array_to_fits(output_norm_SN_rss_fits_filename, norm_SN_rss__yx, overwrite=True) array_to_fits(output_area_rss_fits_filename, area_rss__yx, overwrite=True) seg_cube_filename = f'SEG.cube.{name}.fits' array_to_fits(seg_cube_filename, cube__wyx, header=h_set_cube, overwrite=True) seg_img_filename = f'SEG_img_{name}.fits' array_to_fits(seg_img_filename, V_slice__yx, overwrite=True) print_done(time_ini=time_ini_loc, message='FoV segmentation DONE!') ########################################################################## if config.getboolean('debug', False): return ########################################################################## # Begin RSS analysis # # time_ini_loc = print_block_init('RSS analysis...') # RSS analysis with sigma guided (disp_min) time_ini_auto_ssp = print_block_init('auto_ssp_elines_rnd_rss(): Analyzing RSS with guided sigma...') tmp_cf = input_config['auto_ssp_rss_sigma_guided'] auto_ssp_rss_output_filename = tmp_cf.get('output_filename', 'auto_ssp.CS_few.NAME.rss.out').replace('NAME', name) remove_isfile(auto_ssp_rss_output_filename) auto_ssp_rss_output_elines_filename = 'elines_' + auto_ssp_rss_output_filename remove_isfile(auto_ssp_rss_output_elines_filename) auto_ssp_rss_output_coeffs_filename = 'coeffs_' + auto_ssp_rss_output_filename remove_isfile(auto_ssp_rss_output_coeffs_filename) auto_ssp_output_fits_filename = 'output.' + auto_ssp_rss_output_filename + '.fits.gz' remove_isfile(auto_ssp_output_fits_filename) with open(auto_ssp_rss_output_elines_filename,"w") as elines_out, open(auto_ssp_rss_output_coeffs_filename,"w") as coeffs_out, open(auto_ssp_rss_output_filename,"w") as summary_out: r = auto_ssp_elines_rnd_rss_main( # wavelength=copy(org_wave__w), rss_flux=copy(org_rss__sw), rss_eflux=copy(org_e_rss__sw), wavelength=copy(org_wave__w), rss_flux=org_rss__sw, rss_eflux=org_e_rss__sw, output_files=(elines_out,coeffs_out,summary_out), ssp_file=tmp_cf.get('models_file'), ssp_nl_fit_file=tmp_cf.get('nonlinear_models_file'), sigma_inst=tmp_cf.getfloat('instrumental_dispersion'), out_file=auto_ssp_rss_output_filename, config_file=config_file, mask_list=tmp_cf.get('mask_file', None), elines_mask_file=tmp_cf.get('emission_lines_mask_file', None), min=-2, max=5, w_min=w_min, w_max=w_max, nl_w_min=nl_w_min, nl_w_max=nl_w_max, input_redshift=red, delta_redshift=d_red, min_redshift=min_red, max_redshift=max_red, input_sigma=guess_sigma, delta_sigma=tmp_cf.getfloat('delta_sigma'), min_sigma=tmp_cf.getfloat('min_sigma'), max_sigma=max_sigma, input_AV=tmp_cf.getfloat('input_AV'), delta_AV=tmp_cf.getfloat('delta_AV'), min_AV=tmp_cf.getfloat('min_AV'), max_AV=tmp_cf.getfloat('max_AV'), R_V=tmp_cf.getfloat('RV'), extlaw=tmp_cf.get('extlaw'), plot=plot, is_guided_sigma=True, fit_sigma_rnd=tmp_cf.getboolean('fit_sigma_rnd', True), ) model_spectra, results, results_coeffs, output_el_models = r # write FITS RSS output output_rss__tsw = dump_rss_output(out_file_fit=auto_ssp_output_fits_filename, wavelength=org_wave__w, model_spectra=model_spectra ) print_done(time_ini=time_ini_auto_ssp, message='auto_ssp_elines_rnd_rss(): DONE!') _tmpcf = ConfigAutoSSP(config_file) k = f'{_tmpcf.systems[0]["start_w"]}_{_tmpcf.systems[0]["end_w"]}' index_model_Ha = 0 disp_max_elines = output_el_models[k]['sigma'][index_model_Ha][0]*1.5 del(_tmpcf) # CLEAN MEMORY del(org_rss__sw) del(output_el_models) del(results_coeffs) gc.collect() # # EL: missing plot scripts plot_output_auto_ssp_elines_several_Av_log_rss_all.pl # # creating GAS cube time_init_gas_cube = print_block_init('Creating GAS cube and RSS...') # input_auto_ssp_rss_output_filename = f'output.{auto_ssp_rss_output_filename}.fits.gz' # output_rss__tsw, h_output_rss = get_data_from_fits(input_auto_ssp_rss_output_filename, header=True) print('# -----=> FIT3D_output_rss_seg2cube()', flush=True) ssp_mod_tmp__wyx = rss_seg2cube(copy(output_rss__tsw[1]), seg_map__yx) ssp_mod_cube__wyx = ssp_mod_tmp__wyx*scale_seg__yx tmp_cube__wyx = org_cube__wyx - ssp_mod_cube__wyx # smooth print('# -----=> smooth_spec_clip_cube()', flush=True) smooth_cube__wyx = smooth_spec_clip_cube(copy(tmp_cube__wyx), wavebox_width=75, sigma=1.5, wavepix_min=10, wavepix_max=1860) # generate GAS rss input print('# -----=> spec_extract_cube_mean()', flush=True) gas_cube__wyx = tmp_cube__wyx - smooth_cube__wyx gas_rss__sw, _, _x_mean, _y_mean, npt_mean = spec_extract_cube_mean(copy(gas_cube__wyx), seg_map__yx) # spec_extract_cube_mean(gas_cube_filename, cont_seg_filename, gas_seg_rss_filename) ssp_mod_filename = config_seg.get('ssp_mod_file', 'SSP_mod.NAME.cube.fits').replace('NAME', name) array_to_fits(ssp_mod_filename, ssp_mod_cube__wyx, overwrite=True) gas_cube_filename = config_seg.get('gas_cube_file', 'GAS.NAME.cube.fits').replace('NAME', name) array_to_fits(gas_cube_filename, gas_cube__wyx, header=h_set_cube, overwrite=True) print(f'# {gas_cube_filename} created') gas_seg_rss_filename = config_seg.get('gas_seg_file', 'GAS.CS.NAME.RSS.fits').replace('NAME', name) array_to_fits(gas_seg_rss_filename, gas_rss__sw, header=h_set_rss, overwrite=True) size = np.sqrt(nx**2+ny**2)/(2*ns) output_rss_txt = gas_seg_rss_filename.replace('fits', 'pt.txt') with open(output_rss_txt, 'w') as f: output_header = f'C {size} {size} 0' print(output_header, file=f) # output_header = f'(1) id\n(2) S/N\n(3) Signal\n(4) Noise' np.savetxt(f, list(zip(list(range(ns)), _x_mean/npt_mean, _y_mean/npt_mean, [1]*ns)), fmt=['%d'] + 2*['%.18g'] + ['%d'], delimiter=' ') print(f'# {gas_seg_rss_filename} and {output_rss_txt} created') ssp_mod_tmp_filename = f'SSP_mod_tmp.{name}.cube.fits' array_to_fits(ssp_mod_tmp_filename, ssp_mod_tmp__wyx, header=h_set_cube, overwrite=True) if save_aux_files: # EL: write TMP cube is not needed anymore, but is present as a sanity check tmp_cube_filename = f'TMP.{name}.cube.fits' array_to_fits(tmp_cube_filename, tmp_cube__wyx, overwrite=True) smooth_cube_filename = f'smooth.{name}.cube.fits' array_to_fits(smooth_cube_filename, smooth_cube__wyx, header=h_set_cube, overwrite=True) print_done(time_ini=time_init_gas_cube) del(ssp_mod_tmp__wyx) del(tmp_cube__wyx) del(gas_rss__sw) gc.collect() # analysis of SSP spectra (guided RSS) time_ini_auto_ssp = print_block_init('auto_ssp_elines_rnd_rss(): Analyzing RSS with guided non-linear fit...') tmp_cf = input_config['auto_ssp_rss_guided'] auto_ssp_rss_guided_output_filename = tmp_cf.get('output_filename', 'auto_ssp.CS.NAME.rss.out').replace('NAME', name) remove_isfile(auto_ssp_rss_guided_output_filename) auto_ssp_rss_guided_output_elines_filename = 'elines_' + auto_ssp_rss_guided_output_filename remove_isfile(auto_ssp_rss_guided_output_elines_filename) auto_ssp_rss_guided_output_coeffs_filename = 'coeffs_' + auto_ssp_rss_guided_output_filename remove_isfile(auto_ssp_rss_guided_output_coeffs_filename) auto_ssp_output_fits_filename = 'output.' + auto_ssp_rss_guided_output_filename + '.fits.gz' remove_isfile(auto_ssp_output_fits_filename) guided_nl = True AV = results[5] e_AV = results[6] redshift = results[7] e_redshift = results[8] sigma = results[9] e_sigma = results[10] # AV, e_AV, redshift, e_redshift, sigma, e_sigma = np.genfromtxt( # fname=auto_ssp_rss_output_filename, # comments="#", delimiter=",", usecols=(5,6,7,8,9,10), unpack=True # ) input_guided = (redshift, sigma, AV) input_guided_errors = (e_redshift, e_sigma, e_AV) # SSP spectra ssp_rss__tsw = output_rss__tsw[0] - (output_rss__tsw[2] - output_rss__tsw[1]) with open(auto_ssp_rss_guided_output_elines_filename,"w") as elines_out, open(auto_ssp_rss_guided_output_coeffs_filename,"w") as coeffs_out, open(auto_ssp_rss_guided_output_filename,"w") as summary_out: r = auto_ssp_elines_rnd_rss_main( # wavelength=copy(org_wave__w), rss_flux=copy(ssp_rss__tsw), rss_eflux=copy(org_e_rss__sw), wavelength=copy(org_wave__w), rss_flux=ssp_rss__tsw, rss_eflux=org_e_rss__sw, output_files=(elines_out,coeffs_out,summary_out), input_guided=input_guided, input_guided_errors=input_guided_errors, ssp_file=tmp_cf.get('models_file'), ssp_nl_fit_file=tmp_cf.get('nonlinear_models_file'), sigma_inst=tmp_cf.getfloat('instrumental_dispersion'), out_file=auto_ssp_rss_guided_output_filename, config_file=tmp_cf.get('config_file'), mask_list=tmp_cf.get('mask_file', None), elines_mask_file=tmp_cf.get('emission_lines_mask_file', None), min=-2, max=5, w_min=w_min, w_max=w_max, nl_w_min=nl_w_min, nl_w_max=nl_w_max, input_redshift=red, delta_redshift=d_red, min_redshift=min_red, max_redshift=max_red, input_sigma=guess_sigma, delta_sigma=tmp_cf.getfloat('delta_sigma'), min_sigma=tmp_cf.getfloat('min_sigma'), max_sigma=max_sigma, input_AV=tmp_cf.getfloat('input_AV'), delta_AV=tmp_cf.getfloat('delta_AV'), min_AV=tmp_cf.getfloat('min_AV'), max_AV=tmp_cf.getfloat('max_AV'), R_V=tmp_cf.getfloat('RV'), extlaw=tmp_cf.get('extlaw'), plot=plot, fit_sigma_rnd=tmp_cf.getboolean('fit_sigma_rnd', True), ) del(output_rss__tsw) del(org_e_rss__sw) del(results) gc.collect() model_spectra, results, results_coeffs, _ = r # write FITS RSS output output_guided_rss__tsw = dump_rss_output(out_file_fit=auto_ssp_output_fits_filename, wavelength=org_wave__w, model_spectra=model_spectra ) print_done(time_ini=time_ini_auto_ssp, message='auto_ssp_elines_rnd_rss(): DONE!') print_done(time_ini=time_ini_loc, message='RSS analysis DONE!') # copy elines from auto_ssp disp_min run to final name. shutil.copy( src=f'elines_{auto_ssp_rss_output_filename}', dst=f'elines_{auto_ssp_rss_guided_output_filename}' ) # # EL: missing plot scripts plot_output_auto_ssp_elines_several_Av_log_rss_all.pl # ########################################################################## ########################################################################## # Indices analysis # time_ini_ind_seg = print_block_init('Indices analysis... (CS)') print('# -----=> get_index_output_auto_ssp_elines_rnd_rss_outFIT3D()', flush=True) tmp_cf = input_config['indices'] indices_output_filename = tmp_cf.get('output_file', 'indices.CS.NAME.rss.out').replace('NAME', name) redshift__s = results[_SSP_OUTPUT_INDEX['redshift']] indices, indices__Iyx = get_index( wave__w=org_wave__w, flux_ssp__sw=output_guided_rss__tsw[0] - output_guided_rss__tsw[3], res__sw=output_guided_rss__tsw[4], redshift__s=redshift__s, n_sim=tmp_cf.getint('n_MonteCarlo', 5), plot=plot, seg__yx=seg_map__yx, ) indices_names = list(__indices__.keys()) f = open(indices_output_filename, 'w') for i_s in range(redshift__s.size): for i in range(len(indices_names)): ind_name = indices_names[i] print(f'{ind_name}', end=' ', file=f) EW_now = indices[ind_name]['EW'][i_s] s_EW_now = indices[ind_name]['sigma_EW'][i_s] print(f' {EW_now} {s_EW_now}', end=' ', file=f) med_flux = indices['SN'][i_s] std_flux = indices['e_SN'][i_s] print(f'SN {med_flux} {std_flux}', file=f) f.close() indices_output_fits_filename = tmp_cf.get('output_cube_file', 'indices.CS.NAME.cube.fits').replace('NAME', name) h_indices = { 'COMMENT': 'FIT-header', 'FILENAME': indices_output_fits_filename, } ind_k = [x for x in indices.keys() if 'SN' not in x] n_ind = len(ind_k) for i, k in enumerate(ind_k): h_indices[f'INDEX{i}'] = k h_indices[f'INDEX{i + n_ind + 1}'] = f'e_{k}' i += 1 k = 'SN' h_indices[f'INDEX{i}'] = k h_indices[f'INDEX{i + n_ind + 1}'] = f'e_{k}' array_to_fits(indices_output_fits_filename, indices__Iyx, header=h_indices, overwrite=True) print_done(time_ini=time_ini_ind_seg) ########################################################################## del(output_guided_rss__tsw) del(model_spectra) gc.collect() ########################################################################## # Generating maps # time_ini_loc = print_block_init('Generating maps...') print('# -----=> map_auto_ssp_rnd_seg()', flush=True) input_elines_filename = f'elines_{auto_ssp_rss_guided_output_filename}' map_auto_ssp_rnd_seg( input_elines=input_elines_filename, segmentation_fits=cont_seg_filename, output_prefix=config_seg.get('map_output_prefix', 'map.CS.NAME').replace('NAME', name), inst_disp=config_seg.getfloat('instrumental_dispersion', 0.9), wave_norm=config_seg.getfloat('wave_norm', None), # auto_ssp_config=config_file, ) print_done(time_ini=time_ini_loc) ########################################################################## ########################################################################## # Emission Lines analysis # time_ini_loc = print_block_init('Emission Lines analysis...') wave_ref = config_elines.getfloat('vel_eline_reference_wavelength', __Ha_central_wl__) vel_hr = config_elines.getfloat('vel_eline_half_range', 200) nsearch = config_elines.getint('vel_eline_nsearch', 1) imin = config_elines.getfloat('vel_eline_imin', 0.95) w_min = wave_ref*(1 + min_red - vel_hr/__c__) w_max = wave_ref*(1 + max_red + vel_hr/__c__) print('# -----=> vel_eline_cube()', flush=True) vel_eline_map_filename = config_elines.get('vel_eline_map_file', 've.NAME.vel_map.fits').replace('NAME', name) vel_eline_mask_filename = config_elines.get('vel_eline_mask_map_file', 've.NAME.mask_map.fits').replace('NAME', name) vel_map__yx = np.zeros((ny, nx)) vel_mask_map__yx = np.zeros((ny, nx)) sel_wl = trim_waves(org_wave__w, [w_min, w_max]) wave = org_wave__w[sel_wl] for ixy in itertools.product(range(nx),range(ny)): ix, iy = ixy flux = gas_cube__wyx[sel_wl, iy, ix] ve, mask, npeaks = vel_eline(flux, wave, nsearch=nsearch, imin=imin, wave_ref=wave_ref, set_first_peak=False) vel_map__yx[iy, ix] = ve vel_mask_map__yx[iy, ix] = mask # array_to_fits(vel_eline_map_filename, vel_map__yx, header=h_set_vel, overwrite=True) # h_set_vel['FILENAME'] = vel_eline_mask_filename # array_to_fits(vel_eline_mask_filename, vel_mask_map__yx, header=h_set_vel, overwrite=True) # vel_eline_cube(gas_cube_filename, nsearch=1, imin=0.95, # outfile=vel_eline_prefix, # wave_ref=wave_ref, wmin=w_min, wmax=w_max, # dev='/null') vel_hr = config_elines.getfloat('clean_Ha_map_vel_half_range', 500) min_vel = red*__c__ - vel_hr max_vel = red*__c__ + vel_hr print('# -----=> clean_Ha_map()', flush=True) vel_map_clean__yx, vel_mask_map_clean__yx = clean_Ha_map(copy(vel_map__yx), copy(vel_mask_map__yx), min_vel=min_vel, max_vel=max_vel) h_set_vel = { 'COMMENT': 'vel_eline_cube result', 'WMIN': w_min, 'WMAX': w_max, 'WREF': wave_ref, 'FILENAME': vel_eline_map_filename, } array_to_fits(vel_eline_map_filename, vel_map_clean__yx, overwrite=True, header=h_set_vel) h_set_vel['FILENAME'] = vel_eline_mask_filename array_to_fits(vel_eline_mask_filename, vel_mask_map_clean__yx, overwrite=True, header=h_set_vel) input_config_file = config_elines.get('auto_ssp_config_file') # Preparing config to kin_cube_elines cf = ConfigAutoSSP(input_config_file) vel_hr = config_elines.getfloat('vel_half_range', 300) vel_min = vel - vel_hr vel_max = vel + vel_hr vel_map = copy(vel_map_clean__yx) vel_mask_map = copy(vel_mask_map_clean__yx) sigma_map = None sigma_fixed = None for i_s in range(cf.n_systems): syst = cf.systems[i_s] z_fact = 1 + red w_min = syst['start_w'] w_max = syst['end_w'] w_min_corr = int(w_min*z_fact) w_max_corr = int(w_max*z_fact) elcf = cf.systems_config[i_s] # v0 elcf.guess[0][_MODELS_ELINE_PAR['v0']] = vel elcf.to_fit[0][_MODELS_ELINE_PAR['v0']] = 1 elcf.pars_0[0][_MODELS_ELINE_PAR['v0']] = vel_min elcf.pars_1[0][_MODELS_ELINE_PAR['v0']] = vel_max elcf.links[0][_MODELS_ELINE_PAR['v0']] = -1 # sigma elcf.guess[0][_MODELS_ELINE_PAR['sigma']] = config_elines.getfloat('dispersion_guess', 1.2) elcf.to_fit[0][_MODELS_ELINE_PAR['sigma']] = 1 elcf.pars_0[0][_MODELS_ELINE_PAR['sigma']] = config_elines.getfloat('dispersion_min', 0.5) elcf.pars_1[0][_MODELS_ELINE_PAR['sigma']] = disp_max_elines elcf.links[0][_MODELS_ELINE_PAR['sigma']] = -1 if save_aux_files: # print config to file fname = f'tmp.{name}.{i_s}.config' elcf.print(filename=fname) w_min = int(w_min) w_max = int(w_max) w_minmax_prefix = f'{w_min}_{w_max}' _m_str = 'model' if elcf.n_models > 1: _m_str += 's' print(f'# analyzing {elcf.n_models} {_m_str} in {w_minmax_prefix.replace("_", "-")} wavelength range') output_filename = config_elines.get('output_file', f'KIN.GAS.prefix.NAME.out').replace('NAME', name) output_filename = output_filename.replace('prefix', w_minmax_prefix) f_log = None output_log_filename = config_elines.get('log_file', None) if output_log_filename is not None: output_log_filename = output_log_filename.replace('NAME', name) output_log_filename = output_log_filename.replace('prefix', w_minmax_prefix) f_log = open(output_log_filename, 'w') map_output_prefix = config_elines.get('map_output_prefix', f'map.prefix.NAME').replace('NAME', name).replace('prefix', w_minmax_prefix) map_output_prefix = map_output_prefix.replace('prefix', w_minmax_prefix) fix_disp = 0 if w_min == 3700: sigma_fixed = 1 print('# -----=> kin_cube_elines_rnd()', flush=True) time_ini_kin = print_time(print_seed=False) sel_wl_range = trim_waves(org_wave__w, [w_min_corr, w_max_corr]) if sel_wl_range.astype('int').sum() == 0: print('[kin_cube_elines]: n_wavelength = 0: No avaible spectra to perform the configured analysis.') output_el_models = create_emission_lines_parameters(elcf, shape=(ny, nx)) else: wave_msk__w = org_wave__w[sel_wl_range] flux_msk__wyx = gas_cube__wyx[sel_wl_range] eflux_msk__wyx = None if org_error__wyx is not None: eflux_msk__wyx = copy(org_error__wyx[sel_wl_range]) n_MC = config_elines.getint('n_MonteCarlo', 30) n_loops = config_elines.getint('n_loops', 3) scale_ini = config_elines.getfloat('scale_ini', 0.15) with redirect_stdout(f_log): r = kin_cube_elines_main( wavelength=copy(wave_msk__w), cube_flux=copy(flux_msk__wyx), cube_eflux=eflux_msk__wyx, config=elcf, out_file=output_filename, n_MC=n_MC, n_loops=n_loops, scale_ini=scale_ini, vel_map=vel_map, sigma_map=sigma_map, vel_fixed=0, sigma_fixed=sigma_fixed, redefine_max=1, vel_mask_map=vel_mask_map, memo=0, plot=plot, run_mode=config_elines.get('run_mode', 'RND'), ) if f_log is not None: f_log.close() output_el_spectra, output_el_models = r output_emission_lines_spectra(wave_msk__w, output_el_spectra, h_set_cube, map_output_prefix.replace('map', 'KIN.cube')) output_emission_lines_parameters(map_output_prefix, elcf, output_el_models) print_done(time_ini=time_ini_kin, message='DONE kin_cube_elines_rnd()!') if not i_s: NIIHa_output_el_models = output_el_models print('# -----=> med2df() v_Ha map', flush=True) vel_map_Ha__yx = NIIHa_output_el_models['v0'][index_model_Ha] vel_map_Ha__yx[~np.isfinite(vel_map_Ha__yx)] = 0 vel_map = copy(vel_map_Ha__yx) vel_map_Ha__yx = median_filter(copy(vel_map_Ha__yx), size=(3, 3), mode='reflect') # med2df(f'{map_output_prefix}_vel_00.fits', f'vel.{name}.fits', 3, 3) sigma_map_Ha__yx = NIIHa_output_el_models['sigma'][index_model_Ha] print('# -----=> med2df() disp_Ha map', flush=True) sigma_map_Ha__yx[~np.isfinite(sigma_map_Ha__yx)] = 0 sigma_map = copy(sigma_map_Ha__yx) sigma_map_flux_elines__yx = copy(sigma_map_Ha__yx) disp_map_Ha_filtered__yx = median_filter(sigma_map_Ha__yx*__sigma_to_FWHM__, size=(3, 3), mode='reflect') # med2df(f'{map_output_prefix}_disp_00.fits', f'disp.{name}.fits', 3, 3) print('# -----=> imarith()', flush=True) # imarith(f'disp.{name}.fits', '*', str(__FWHM_to_sigma__), f'disp.{name}.fits') if save_aux_files: vel_filename = f'vel.{name}.fits' remove_isfile(vel_filename) array_to_fits(vel_filename, vel_map_Ha__yx, overwrite=True) disp_filename = f'disp.{name}.fits' remove_isfile(disp_filename) sigma_map_Ha_filtered__yx = __FWHM_to_sigma__*disp_map_Ha_filtered__yx array_to_fits(disp_filename, sigma_map_Ha_filtered__yx, overwrite=True) print(f'# emission lines in wavelength range {w_minmax_prefix.replace("_", "-")} analyzed') print_done(time_ini=time_ini_loc) ########################################################################## ########################################################################## # Momentum analysis of emission lines # time_ini_loc = print_block_init('Momentum analysis of emission lines (flux_elines)...') tmp_cf = input_config['flux_elines'] elines_list = tmp_cf.get('elines_list') flux_elines_output_filename = tmp_cf.get('output_file', 'flux_elines.NAME.cube.fits.gz').replace('NAME', name) print('# -----=> flux_elines_cube_EW()', flush=True) fe_output, fe_header = flux_elines_cube_EW( flux__wyx = gas_cube__wyx, input_header = h_set_cube, n_MC=tmp_cf.getint('n_MonteCarlo', 10), elines_list=elines_list, vel__yx=vel_map_clean__yx, sigma__yx=sigma_map_flux_elines__yx, eflux__wyx=org_error__wyx, flux_ssp__wyx=ssp_mod_cube__wyx, ) array_to_fits(flux_elines_output_filename, fe_output, header=fe_header, overwrite=True) del(fe_output) gc.collect() print_done(time_ini=time_ini_loc) elines_long_list = tmp_cf.get('elines_long_list', None) if elines_long_list is not None: time_ini_loc = print_block_init('Momentum analysis of emission lines (flux_elines_long)...') flux_elines_long_output_filename = flux_elines_output_filename.replace('flux_elines', 'flux_elines_long') print('# -----=> flux_elines_cube_EW()', flush=True) fe_output, fe_header = flux_elines_cube_EW( flux__wyx = gas_cube__wyx, input_header = h_set_cube, n_MC=tmp_cf.getint('n_MonteCarlo', 10), elines_list=elines_long_list, vel__yx=vel_map_clean__yx, sigma__yx=sigma_map_flux_elines__yx, eflux__wyx=org_error__wyx, flux_ssp__wyx=ssp_mod_cube__wyx, ) array_to_fits(flux_elines_long_output_filename, fe_output, header=fe_header, overwrite=True) del(fe_output) gc.collect() print_done(time_ini=time_ini_loc) ########################################################################## del(gas_cube__wyx) gc.collect() ########################################################################## # gzip fits # time_ini_loc = print_block_init('Compressing FITS files...') compresslevel = config.getint('gzip_compresslevel', 6) for file in glob.glob(f'*{name}*.fits'): with open(file, 'rb') as f_in: # EL: compresslevel 6 is biased towars high compression # at expense of speed, also 6 is also the value of gzip # in shell mode. with gzip.open(f'{file}.gz', mode='wb', compresslevel=compresslevel) as f_comp: shutil.copyfileobj(f_in, f_comp) # remove uncompressed file remove(file) print_done(time_ini=time_ini_loc) ########################################################################## ########################################################################## # generate mass age met maps # time_ini_loc = print_block_init('Creating mass, age and metallicities map...') sum_mass_age_filename = config_pack.get('sum_mass_age_out_file_1').replace('NAME', name) print('# -----=> sum_mass_age()', flush=True) sum_mass_age(name, output_csv_filename=sum_mass_age_filename) print_done(time_ini=time_ini_loc) ########################################################################## ########################################################################## ########################################################################## ########################################################################## ########################################################################## ########################################################################## ########################################################################## ########################################################################## # CS binning - SFH fill # The analysis of the datacube spaxel by spaxel deseg_cf = input_config['desegmentation'] deseg_prefix = deseg_cf.get('prefix') w_min, w_max = org_w_minmax nl_w_min, nl_w_max = org_nl_w_minmax ########################################################################## # Begin all spaxels analysis # time_ini_allspx_ana = print_block_init('Begin all spaxels analysis (SFH CS fill)...') auto_ssp_deseg_output_filename = auto_ssp_rss_guided_output_filename.replace('CS', deseg_prefix).replace('rss', 'cube') remove_isfile(auto_ssp_deseg_output_filename) auto_ssp_deseg_output_coeffs_filename = 'coeffs_' + auto_ssp_deseg_output_filename remove_isfile(auto_ssp_deseg_output_coeffs_filename) auto_ssp_deseg_output_coeffs_nhlib_filename = 'coeffs_nhlib_' + auto_ssp_deseg_output_filename remove_isfile(auto_ssp_deseg_output_coeffs_nhlib_filename) auto_ssp_deseg_output_elines_filename = 'elines_' + auto_ssp_deseg_output_filename remove_isfile(auto_ssp_deseg_output_elines_filename) auto_ssp_deseg_output_fits_filename = 'output.' + auto_ssp_deseg_output_filename + '.fits.gz' remove_isfile(auto_ssp_deseg_output_fits_filename) spx_nh_coeffs_filename = 'nh_coeffs.' + auto_ssp_deseg_output_filename + '.fits.gz' remove_isfile(spx_nh_coeffs_filename) with open(auto_ssp_deseg_output_elines_filename,"w") as elines_out, open(auto_ssp_deseg_output_coeffs_filename,"w") as coeffs_out, open(auto_ssp_deseg_output_coeffs_nhlib_filename,"w") as coeffs_nhlib_out, open(auto_ssp_deseg_output_filename,"w") as summary_out: r = auto_ssp_desegmentation_main( # copy(org_wave__w), copy(org_cube__wyx), copy(org_error__wyx), copy(org_wave__w), org_cube__wyx, org_error__wyx, output_files=(elines_out, coeffs_out, coeffs_nhlib_out, summary_out), ssp_file=deseg_cf.get('models_file'), config_file=config_file, out_file=auto_ssp_deseg_output_filename, seg_map=seg_map__yx, rss_nl_pars=(results[7], results[9], results[5]), rss_coeffs=(results_coeffs[0], results_coeffs[1], results_coeffs[3]), ssp_nl_fit_file=deseg_cf.get('nonlinear_models_file'), sigma_inst=deseg_cf.getfloat('instrumental_dispersion'), mask_list=deseg_cf.get('mask_file', None), elines_mask_file=deseg_cf.get('emission_lines_mask_file', None), w_min=w_min, w_max=w_max, nl_w_min=nl_w_min, nl_w_max=nl_w_max, R_V=deseg_cf.getfloat('RV'), extlaw=deseg_cf.get('extlaw'), fit_sigma_rnd=deseg_cf.getboolean('fit_sigma_rnd', True), plot=None, sps_class=StPopSynt, kdtree_nodes=deseg_cf.getint('kdtree_nodes', 4), # store_models=True, store_results=True, store_el=True, store_coeffs=True, ) del(results) del(results_coeffs) gc.collect() model_spectra, results, results_coeffs, results_coeffs_nhlib, spx_nh_coeffs, output_el_models, analisys_sequence = r analisys_sequence_filename = f'{deseg_prefix}.{name}.analisys_sequence.fits.gz' array_to_fits(analisys_sequence_filename, np.asarray(analisys_sequence), overwrite=True) dump_rss_output(out_file_fit=auto_ssp_deseg_output_fits_filename, wavelength=org_wave__w, model_spectra=model_spectra) # dump_auto_ssp_el_output(output_el_models, config_file, deseg_prefix) array_to_fits(spx_nh_coeffs_filename, spx_nh_coeffs, overwrite=True) print_done(time_ini=time_ini_allspx_ana, message=f'{name}: all spaxels analysis (SFH CS fill): DONE!') _tmpcf = ConfigAutoSSP(config_file) k = f'{_tmpcf.systems[0]["start_w"]}_{_tmpcf.systems[0]["end_w"]}' index_model_Ha = 0 _ix0, _iy0 = analisys_sequence[0] disp_max_elines = output_el_models[k]['sigma'][index_model_Ha][_iy0, _ix0]*1.5 del(_tmpcf) del(output_el_models) map_auto_ssp_deseg( name=name, wavelength=org_wave__w, model_spectra=model_spectra, results=results, results_coeffs=results_coeffs, ssp_file=input_config['auto_ssp_rss_guided'].get('models_file'), V__yx=mV__yx, mask_map__yx=mask_map__yx, prefix=deseg_prefix, output_prefix=config_seg.get('map_output_prefix').replace('CS', deseg_prefix).replace('NAME', name), wave_norm=config_seg.getfloat('wave_norm', None), inst_disp=config_seg.getfloat('instrumental_dispersion'), sum_mass_age_out_file=f'{deseg_prefix}.' + config_pack.get('sum_mass_age_out_file_1', f'MASS.NAME.csv').replace('NAME', name), cont_seg_file=f'{deseg_prefix}.' + config_seg.get('cont_seg_file').replace('NAME', name) + '.gz', scale_seg_file=f'{deseg_prefix}.' + config_seg.get('scale_seg_file').replace('NAME', name) + '.gz', ) ########################################################################## # SFH CS fill Indices analysis # time_ini_ind_xy = print_block_init('Indices analysis... (XY)') print('# -----=> get_index_output_auto_ssp_elines_rnd_rss_outFIT3D()', flush=True) indices_deseg_output_filename = indices_output_filename.replace('CS', deseg_prefix).replace('rss', 'cube') redshift__yx = results[7] indices, indices__Iyx = get_index_cube( wave__w=org_wave__w, flux_ssp__wyx=model_spectra[0] - model_spectra[3], res__wyx=model_spectra[4], redshift__yx=redshift__yx, n_sim=tmp_cf.getint('n_MonteCarlo', 5), plot=plot, ) indices_names = list(__indices__.keys()) f = open(indices_deseg_output_filename, 'w') for ixy in itertools.product(range(nx), range(ny)): i_x, i_y = ixy for i in range(len(indices_names)): ind_name = indices_names[i] print(f'{ind_name}', end=' ', file=f) EW_now = indices[ind_name]['EW'][i_y, i_x] s_EW_now = indices[ind_name]['sigma_EW'][i_y, i_x] print(f' {EW_now} {s_EW_now}', end=' ', file=f) med_flux = indices['SN'][i_y, i_x] std_flux = indices['e_SN'][i_y, i_x] print(f'SN {med_flux} {std_flux}', file=f) f.close() indices_deseg_output_fits_filename = indices_output_fits_filename.replace('CS', deseg_prefix) h_indices = { 'COMMENT': 'FIT-header', 'FILENAME': indices_deseg_output_fits_filename, } ind_k = [x for x in indices.keys() if 'SN' not in x] n_ind = len(ind_k) for i, k in enumerate(ind_k): h_indices[f'INDEX{i}'] = k h_indices[f'INDEX{i + n_ind + 1}'] = f'e_{k}' i += 1 k = 'SN' h_indices[f'INDEX{i}'] = k h_indices[f'INDEX{i + n_ind + 1}'] = f'e_{k}' array_to_fits(indices_deseg_output_fits_filename, indices__Iyx, header=h_indices, overwrite=True) print_done(time_ini=time_ini_ind_xy) ########################################################################## ########################################################################## # Creating gas cube ########################################################################## ssp_mod_cube__wyx = model_spectra[1] gas_cube__wyx = org_cube__wyx - ssp_mod_cube__wyx gas_cube_filename = gas_cube_filename.replace(name, f'{deseg_prefix}.{name}') array_to_fits(gas_cube_filename, gas_cube__wyx, header=h_set_cube, overwrite=True) print(f'# {gas_cube_filename} created') ########################################################################## # Emission Lines analysis # time_ini_loc = print_block_init('Emission Lines analysis...') wave_ref = config_elines.getfloat('vel_eline_reference_wavelength', __Ha_central_wl__) vel_hr = config_elines.getfloat('vel_eline_half_range', 200) nsearch = config_elines.getint('vel_eline_nsearch', 1) imin = config_elines.getfloat('vel_eline_imin', 0.95) w_min = wave_ref*(1 + min_red - vel_hr/__c__) w_max = wave_ref*(1 + max_red + vel_hr/__c__) print('# -----=> vel_eline_cube()', flush=True) vel_eline_map_filename = vel_eline_map_filename.replace(name, f'{deseg_prefix}.{name}') vel_eline_mask_filename = vel_eline_mask_filename.replace(name, f'{deseg_prefix}.{name}') vel_map__yx = np.zeros((ny, nx)) vel_mask_map__yx = np.zeros((ny, nx)) sel_wl = trim_waves(org_wave__w, [w_min, w_max]) wave = org_wave__w[sel_wl] for ixy in itertools.product(range(nx), range(ny)): ix, iy = ixy flux = gas_cube__wyx[sel_wl, iy, ix] ve, mask, npeaks = vel_eline(flux, wave, nsearch=nsearch, imin=imin, wave_ref=wave_ref, set_first_peak=False) vel_map__yx[iy, ix] = ve vel_mask_map__yx[iy, ix] = mask vel_hr = config_elines.getfloat('clean_Ha_map_vel_half_range', 500) min_vel = red*__c__ - vel_hr max_vel = red*__c__ + vel_hr print('# -----=> clean_Ha_map()', flush=True) vel_map_clean__yx, vel_mask_map_clean__yx = clean_Ha_map(copy(vel_map__yx), copy(vel_mask_map__yx), min_vel=min_vel, max_vel=max_vel) h_set_vel = { 'COMMENT': 'vel_eline_cube result', 'WMIN': w_min, 'WMAX': w_max, 'WREF': wave_ref, 'FILENAME': vel_eline_map_filename, } array_to_fits(vel_eline_map_filename, vel_map_clean__yx, overwrite=True, header=h_set_vel) h_set_vel['FILENAME'] = vel_eline_mask_filename array_to_fits(vel_eline_mask_filename, vel_mask_map_clean__yx, overwrite=True, header=h_set_vel) input_config_file = config_elines.get('auto_ssp_config_file') # Preparing config to kin_cube_elines cf = ConfigAutoSSP(input_config_file) vel_hr = config_elines.getfloat('vel_half_range', 300) vel_min = vel - vel_hr vel_max = vel + vel_hr vel_map = copy(vel_map_clean__yx) vel_mask_map = copy(vel_mask_map_clean__yx) sigma_map = None sigma_fixed = None for i_s in range(cf.n_systems): syst = cf.systems[i_s] z_fact = 1 + red w_min = syst['start_w'] w_max = syst['end_w'] w_min_corr = int(w_min*z_fact) w_max_corr = int(w_max*z_fact) elcf = cf.systems_config[i_s] # v0 elcf.guess[0][_MODELS_ELINE_PAR['v0']] = vel elcf.to_fit[0][_MODELS_ELINE_PAR['v0']] = 1 elcf.pars_0[0][_MODELS_ELINE_PAR['v0']] = vel_min elcf.pars_1[0][_MODELS_ELINE_PAR['v0']] = vel_max elcf.links[0][_MODELS_ELINE_PAR['v0']] = -1 # sigma elcf.guess[0][_MODELS_ELINE_PAR['sigma']] = config_elines.getfloat('dispersion_guess', 1.2) elcf.to_fit[0][_MODELS_ELINE_PAR['sigma']] = 1 elcf.pars_0[0][_MODELS_ELINE_PAR['sigma']] = config_elines.getfloat('dispersion_min', 0.5) elcf.pars_1[0][_MODELS_ELINE_PAR['sigma']] = disp_max_elines elcf.links[0][_MODELS_ELINE_PAR['sigma']] = -1 if save_aux_files: # print config to file fname = f'tmp.{deseg_prefix}.{name}.{i_s}.config' elcf.print(filename=fname) w_min = int(w_min) w_max = int(w_max) w_minmax_prefix = f'{w_min}_{w_max}' _m_str = 'model' if elcf.n_models > 1: _m_str += 's' print(f'# analyzing {elcf.n_models} {_m_str} in {w_minmax_prefix.replace("_", "-")} wavelength range') output_filename = config_elines.get('output_file', f'KIN.GAS.prefix.NAME.out').replace('NAME', f'{deseg_prefix}.{name}') output_filename = output_filename.replace('prefix', w_minmax_prefix) f_log = None output_log_filename = config_elines.get('log_file', None) if output_log_filename is not None: output_log_filename = output_log_filename.replace('NAME', name) output_log_filename = output_log_filename.replace('prefix', w_minmax_prefix) f_log = open(output_log_filename, 'w') map_output_prefix = config_elines.get('map_output_prefix', f'map.prefix.NAME').replace('NAME', f'{deseg_prefix}.{name}').replace('prefix', w_minmax_prefix) map_output_prefix = map_output_prefix.replace('prefix', w_minmax_prefix) fix_disp = 0 if w_min == 3700: sigma_fixed = 1 print('# -----=> kin_cube_elines_rnd()', flush=True) time_ini_kin = print_time(print_seed=False) sel_wl_range = trim_waves(org_wave__w, [w_min_corr, w_max_corr]) if sel_wl_range.astype('int').sum() == 0: print('[kin_cube_elines]: n_wavelength = 0: No avaible spectra to perform the configured analysis.') output_el_models = create_emission_lines_parameters(elcf, shape=(ny, nx)) else: wave_msk__w = org_wave__w[sel_wl_range] flux_msk__wyx = gas_cube__wyx[sel_wl_range] eflux_msk__wyx = None if org_error__wyx is not None: eflux_msk__wyx = copy(org_error__wyx[sel_wl_range]) n_MC = config_elines.getint('n_MonteCarlo', 30) n_loops = config_elines.getint('n_loops', 3) scale_ini = config_elines.getfloat('scale_ini', 0.15) with redirect_stdout(f_log): r = kin_cube_elines_main( wavelength=copy(wave_msk__w), cube_flux=copy(flux_msk__wyx), cube_eflux=eflux_msk__wyx, config=elcf, out_file=output_filename, n_MC=n_MC, n_loops=n_loops, scale_ini=scale_ini, vel_map=vel_map, sigma_map=sigma_map, vel_fixed=0, sigma_fixed=sigma_fixed, redefine_max=1, vel_mask_map=vel_mask_map, memo=0, plot=plot, run_mode=config_elines.get('run_mode', 'RND'), ) if f_log is not None: f_log.close() output_el_spectra, output_el_models = r output_emission_lines_spectra(wave_msk__w, output_el_spectra, h_set_cube, map_output_prefix.replace('map', 'KIN.cube')) output_emission_lines_parameters(map_output_prefix, elcf, output_el_models) print_done(time_ini=time_ini_kin, message='DONE kin_cube_elines_rnd()!') if not i_s: NIIHa_output_el_models = output_el_models print('# -----=> med2df() v_Ha map', flush=True) vel_map_Ha__yx = NIIHa_output_el_models['v0'][index_model_Ha] vel_map_Ha__yx[~np.isfinite(vel_map_Ha__yx)] = 0 vel_map = copy(vel_map_Ha__yx) vel_map_Ha__yx = median_filter(copy(vel_map_Ha__yx), size=(3, 3), mode='reflect') # med2df(f'{map_output_prefix}_vel_00.fits', f'vel.{name}.fits', 3, 3) sigma_map_Ha__yx = NIIHa_output_el_models['sigma'][index_model_Ha] print('# -----=> med2df() disp_Ha map', flush=True) sigma_map_Ha__yx[~np.isfinite(sigma_map_Ha__yx)] = 0 sigma_map = copy(sigma_map_Ha__yx) sigma_map_flux_elines__yx = copy(sigma_map_Ha__yx) disp_map_Ha_filtered__yx = median_filter(sigma_map_Ha__yx*__sigma_to_FWHM__, size=(3, 3), mode='reflect') # med2df(f'{map_output_prefix}_disp_00.fits', f'disp.{name}.fits', 3, 3) print('# -----=> imarith()', flush=True) # imarith(f'disp.{name}.fits', '*', str(__FWHM_to_sigma__), f'disp.{name}.fits') if save_aux_files: vel_filename = f'vel.{deseg_prefix}.{name}.fits' remove_isfile(vel_filename) array_to_fits(vel_filename, vel_map_Ha__yx, overwrite=True) disp_filename = f'disp.{deseg_prefix}.{name}.fits' remove_isfile(disp_filename) sigma_map_Ha_filtered__yx = __FWHM_to_sigma__*disp_map_Ha_filtered__yx array_to_fits(disp_filename, sigma_map_Ha_filtered__yx, overwrite=True) print(f'# emission lines in wavelength range {w_minmax_prefix.replace("_", "-")} analyzed') print_done(time_ini=time_ini_loc) ########################################################################## ########################################################################## # Momentum analysis of emission lines # time_ini_loc = print_block_init('Momentum analysis of emission lines (flux_elines)...') tmp_cf = input_config['flux_elines'] elines_list = tmp_cf.get('elines_list') flux_elines_deseg_output_filename = tmp_cf.get('output_file', 'flux_elines.NAME.cube.fits.gz').replace('NAME', f'{deseg_prefix}.{name}') print('# -----=> flux_elines_cube_EW()', flush=True) fe_output, fe_header = flux_elines_cube_EW( flux__wyx = gas_cube__wyx, input_header = h_set_cube, n_MC=tmp_cf.getint('n_MonteCarlo', 10), elines_list=elines_list, vel__yx=vel_map_clean__yx, sigma__yx=sigma_map_flux_elines__yx, eflux__wyx=org_error__wyx, flux_ssp__wyx=ssp_mod_cube__wyx, ) array_to_fits(flux_elines_deseg_output_filename, fe_output, header=fe_header, overwrite=True) del(fe_output) gc.collect() print_done(time_ini=time_ini_loc) elines_long_list = tmp_cf.get('elines_long_list', None) if elines_long_list is not None: time_ini_loc = print_block_init('Momentum analysis of emission lines (flux_elines_long)...') flux_elines_long_deseg_output_filename = flux_elines_deseg_output_filename.replace('flux_elines', 'flux_elines_long') print('# -----=> flux_elines_cube_EW()', flush=True) fe_output, fe_header = flux_elines_cube_EW( flux__wyx = gas_cube__wyx, input_header = h_set_cube, n_MC=tmp_cf.getint('n_MonteCarlo', 10), elines_list=elines_long_list, vel__yx=vel_map_clean__yx, sigma__yx=sigma_map_flux_elines__yx, eflux__wyx=org_error__wyx, flux_ssp__wyx=ssp_mod_cube__wyx, ) array_to_fits(flux_elines_long_deseg_output_filename, fe_output, header=fe_header, overwrite=True) del(fe_output) gc.collect() print_done(time_ini=time_ini_loc) ########################################################################## # gzip fits # time_ini_loc = print_block_init('Compressing FITS files...') compresslevel = config.getint('gzip_compresslevel', 6) for file in glob.glob(f'*{name}*.fits'): with open(file, 'rb') as f_in: # EL: compresslevel 6 is biased towars high compression # at expense of speed, also 6 is also the value of gzip # in shell mode. with gzip.open(f'{file}.gz', mode='wb', compresslevel=compresslevel) as f_comp: shutil.copyfileobj(f_in, f_comp) # remove uncompressed file remove(file) print_done(time_ini=time_ini_loc) ########################################################################## ########################################################################## ########################################################################## ########################################################################## ########################################################################## ########################################################################## ########################################################################## ########################################################################## ########################################################################## ########################################################################## # Creating output datacubes ########################################################################## ########################################################################## ########################################################################## # Final pack, and store files # # add a / at the end of directories names time_ini_final_pack = print_block_init('Begin final pack of files...') if not dir_out.endswith('/'): dir_out += '/' if not dir_datacubes_out.endswith('/'): dir_datacubes_out += '/' ########################################################################## ########################################################################## # Mask maps # time_ini_loc = print_block_init('Masking maps...') for fname in glob.glob(f'map.CS*{name}*.fits.gz'): print(f' ---=> imarith() {fname}', flush=True) imarith(fname, '*', f'{SN_mask_map_filename}.gz', fname) for fname in glob.glob(f'map.{deseg_prefix}*{name}*.fits.gz'): print(f' ---=> imarith() {fname}', flush=True) imarith(fname, '*', f'{SN_mask_map_filename}.gz', fname) print_done(time_ini=time_ini_loc) ########################################################################## ########################################################################## # Pack files # time_ini_loc = print_block_init('Packing files...') print('# -----=> pack_NAME()', flush=True) pack_NAME( name, ssp_pack_filename=config_pack.get('SSP_file'), elines_pack_filename=config_pack.get('ELINES_file'), sfh_pack_filename=config_pack.get('SFH_file'), mass_out_filename=config_pack.get('sum_mass_age_out_file_2').replace('NAME', name) ) # # SFH CS fill datacubes # new_SFH_pack, new_SSP_pack, new_ELINES_pack = create_deseg_packs( name=name, org_SFH_pack=config_pack.get('SFH_file'), org_SSP_pack=config_pack.get('SSP_file'), org_ELINES_pack=config_pack.get('ELINES_file'), org_output_prefix=config_seg.get('map_output_prefix'), org_cont_seg=config_seg.get('cont_seg_file'), org_scale_seg=config_seg.get('scale_seg_file'), prefix=deseg_prefix, ) pack_results_name(name, new_SFH_pack, f'SFH.{deseg_prefix}') pack_results_name(name, new_SSP_pack, f'SSP.{deseg_prefix}') pack_results_name(name, new_ELINES_pack, f'ELINES.{deseg_prefix}') # # pack dataproducts # datacubes_names_str = config_pack.get('datacubes_names') datacubes_mask_name = config_pack.get('mask_name', 'MASK') datacubes_mask_filename = config_pack.get('mask_file', None) if datacubes_mask_filename is not None: datacubes_mask_filename = datacubes_mask_filename.replace('NAME', name) datacubes_names_str = datacubes_names_str.replace(',\'MASK\',', f',\'{datacubes_mask_name}\',') datacubes_names = list(eval(datacubes_names_str)) datacubes_output_filename = config_pack.get('datacubes_output_filename').replace('NAME', name) pack_dp = { # CS 'SSP': f'{name}.SSP.cube.fits.gz', 'SFH': f'{name}.SFH.cube.fits.gz', 'ELINES': f'{name}.ELINES.cube.fits.gz', 'INDICES': indices_output_fits_filename + '.gz', 'FLUX_ELINES': flux_elines_output_filename, 'FLUX_ELINES_LONG': None if elines_long_list is None else flux_elines_long_output_filename, # deseg 'DESEG_SSP': f'{name}.SSP.{deseg_prefix}.cube.fits.gz', 'DESEG_SFH': f'{name}.SFH.{deseg_prefix}.cube.fits.gz', 'DESEG_ELINES': f'{name}.ELINES.{deseg_prefix}.cube.fits.gz', 'DESEG_INDICES': indices_deseg_output_fits_filename + '.gz', 'DESEG_FLUX_ELINES': flux_elines_deseg_output_filename, 'DESEG_FLUX_ELINES_LONG': None if elines_long_list is None else flux_elines_long_deseg_output_filename, # masks datacubes_mask_name: datacubes_mask_filename, 'SELECT_REG': SN_mask_map_filename + '.gz', } pack_dataproducts(header=org_h, files_dict=pack_dp, datacubes_names=datacubes_names, output_filename=datacubes_output_filename) print_done(time_ini=time_ini_loc) ########################################################################## ########################################################################## # Create final directory # time_ini_loc = print_block_init('Store generated files...') if 'manga' in name: try: _, plate, ifu = name.split('-') dir_out += f'{plate}/{ifu}/' except: dir_out += f'{name}/' else: dir_out += f'{name}/' try: makedirs(dir_out) except FileExistsError: print(f'{basename(sys.argv[0])}: {dir_out}: directory already exists') try: makedirs(dir_datacubes_out) except FileExistsError: print(f'{basename(sys.argv[0])}: {dir_datacubes_out}: directory already exists') # copy and move files # the filename should be in dst in order to rewrite file print(f'Moving files to {dir_out}', flush=True) for file in glob.glob(f'*{name}*'): shutil.move(src=file, dst=f'{dir_out}{file}') print_done(time_ini=time_ini_loc, message='DONE move files!') print(f'Copying datacubes to {dir_datacubes_out}', flush=True) cen_auto_ssp_output_fits_filename += '.gz' files_to_copy = [ int_auto_ssp_output_filename, cen_auto_ssp_output_fits_filename, int_spec_filename, # CS pack_dp['SSP'], pack_dp['SFH'], pack_dp['ELINES'], pack_dp['FLUX_ELINES'], pack_dp['INDICES'], # deseg pack_dp['DESEG_SSP'], pack_dp['DESEG_SFH'], pack_dp['DESEG_ELINES'], pack_dp['DESEG_FLUX_ELINES'], pack_dp['DESEG_INDICES'], # final datacube datacubes_output_filename, ] if elines_long_list is not None: files_to_copy.append(pack_dp['FLUX_ELINES_LONG']) files_to_copy.append(pack_dp['DESEG_FLUX_ELINES_LONG']) for file in files_to_copy: shutil.copy(src=f'{dir_out}{file}', dst=f'{dir_datacubes_out}') print_done(time_ini=time_ini_loc, message='DONE: store files!') ########################################################################## print_done(time_ini=time_ini_final_pack, message='DONE: final pack!') ########################################################################## # Plot needed # TODO ana_single_MaNGA_plot.pl ########################################################################## return dir_out def main(args): name = args.name force_z = args.redshift x0 = args.x0 y0 = args.y0 # only_binned_spaxels = args.only_binned_spaxels config = configparser.ConfigParser( # allow a variables be set without value allow_no_value=True, # allows duplicated keys in different sections strict=False, # deals with variables inside configuratio file interpolation=configparser.ExtendedInterpolation()) config.read(args.config_file) dir_log = 'log' log_file = config['general'].get('log_file', None) f = None if (log_file is not None): if 'NAME' in log_file: log_file = log_file.replace('NAME', name) omode = 'w' else: omode='a' if log_file[0:5] != '/dev/': try: makedirs(dir_log) log_file = f'{dir_log}/{log_file}' except FileExistsError: print(f'{basename(sys.argv[0])}: {dir_log}: directory already exists') log_file = f'{dir_log}/{log_file}' except: print(f'{basename(sys.argv[0])}: {dir_log}: cannot create directory') else: omode = 'w' f = open(log_file, omode) with redirect_stdout(f): with redirect_stderr(f): init_message = f'Initiating run for galaxy {name} - May the force be with you!' init_message += f'\n# Configuration file = {args.config_file}' init_message += f'\n# Output log file = {log_file}' cube_filename = config['general']['cube_filename'].replace('NAME', name) init_message += f'\n# Cube file = {cube_filename}' time_ini_run = print_block_init(init_message, print_seed=True) seed = config['general'].getint('seed', None) if seed is not None: print(f'# Forcing seed: {seed}') else: seed = time_ini_run np.random.seed(seed) dir_out = ana_single(name, config, force_z=force_z, x0=x0, y0=y0) #, only_binned_spaxels=only_binned_spaxels) print_done(time_ini=time_ini_run, message=f'DONE: galaxy {name} analyzed.') # After some problems with the write in the log_file # we move the log_file after everithing is done. # If the program was not able to create the dir_log directory # in some machines the final writting of the log file # could be incomplete. if (log_file is not None) and os.path.isfile(log_file): f.close() if name in log_file: _tmp = config['general'].get('log_file', 'ana_single.NAME.log').replace('NAME', name) shutil.move(src=log_file, dst=f'{dir_out}{_tmp}') if __name__ == '__main__': main(ReadArgumentsLocal())