From 5eed46afe840c26feac76a0b1ff98f94e8f2155e Mon Sep 17 00:00:00 2001 From: chloe Date: Thu, 24 Mar 2022 14:07:35 -0400 Subject: [PATCH] Read fits from older versions of alf and normalize ext_model --- scripts/read_alf.py | 309 +++++++++++++++++++++++++++++++++++++------- 1 file changed, 260 insertions(+), 49 deletions(-) diff --git a/scripts/read_alf.py b/scripts/read_alf.py index 5841cc2..94d3c15 100644 --- a/scripts/read_alf.py +++ b/scripts/read_alf.py @@ -12,15 +12,24 @@ class Alf(object): def __init__(self, outfiles, read_mcmc=True, info=None, index=False): self.outfiles = outfiles + ## ND: New file reading code, don't need this after rehaul. Delete. #self.legend = info['label'] #self.imf_type = info['imf_type'] - self.nsample = None + #self.nsample = None self.spectra = None + self.model_info = {} if read_mcmc: self.mcmc = np.loadtxt('{0}.mcmc'.format(self.outfiles)) results = ascii.read('{0}.sum'.format(self.outfiles)) + # To keep track of labels for backwards compatibility of + # previous versions of alf + self.filetype = len(results.colnames) + + + ## ND: New file reading code, don't need this after rehaul. Delete. + """ with open('{0}.sum'.format(self.outfiles)) as f: for line in f: if line[0] == '#': @@ -30,18 +39,57 @@ def __init__(self, outfiles, read_mcmc=True, info=None, index=False): self.nchain = float(line.split('=')[1].strip()) elif 'Nsample' in line: self.nsample = float(line.split('=')[1].strip()) + """ - self.labels = np.array([ - 'chi2', 'velz', 'sigma', 'logage', 'zH', 'FeH', 'a', - 'C', 'N', 'Na', 'Mg', 'Si', 'K', 'Ca', 'Ti', 'V', 'Cr', - 'Mn', 'Co', 'Ni', 'Cu', 'Sr', 'Ba', 'Eu', 'Teff', - 'IMF1', 'IMF2', 'logfy', 'sigma2', 'velz2', 'logm7g', - 'hotteff', 'loghot', 'fy_logage', 'logemline_h', - 'logemline_oii', 'logemline_oiii', 'logemline_sii', - 'logemline_ni', 'logemline_nii', 'logtrans', 'jitter', - 'logsky', 'IMF3', 'IMF4', 'h3', 'h4', - 'ML_v','ML_i','ML_k','MW_v', 'MW_i','MW_k' - ]) + with open('{0}.sum'.format(self.outfiles)) as f: + for line in f: + if line[0] == '#': + attr = line.replace('#', '').split('=') + if len(attr) == 2: + try: + self.model_info[attr[0].strip()] = int(attr[1].strip()) + except: + self.model_info[attr[0].strip()] = attr[1].strip() + else: + break + + # Older versions of alf have different labels/different label positions + # Current version has 53 labels + # Older versions have 52 or 50 (i.e. van Dokkum+17 has 50) + if self.filetype == 53: + self.labels = np.array([ + 'chi2', 'velz', 'sigma', 'logage', 'zH', 'FeH', 'a', + 'C', 'N', 'Na', 'Mg', 'Si', 'K', 'Ca', 'Ti', 'V', 'Cr', + 'Mn', 'Co', 'Ni', 'Cu', 'Sr', 'Ba', 'Eu', 'Teff', + 'IMF1', 'IMF2', 'logfy', 'sigma2', 'velz2', 'logm7g', + 'hotteff', 'loghot', 'fy_logage', 'logemline_h', + 'logemline_oii', 'logemline_oiii', 'logemline_sii', + 'logemline_ni', 'logemline_nii', 'logtrans', 'jitter', + 'logsky', 'IMF3', 'IMF4', 'h3', 'h4', + 'ML_v','ML_i','ML_k','MW_v', 'MW_i','MW_k' + ]) + elif self.filetype == 52: + self.labels = np.array([ + 'chi2', 'velz', 'sigma', 'logage', 'zH', 'FeH', 'a', + 'C', 'N', 'Na', 'Mg', 'Si', 'K', 'Ca', 'Ti', 'V', 'Cr', + 'Mn', 'Co', 'Ni', 'Cu', 'Sr', 'Ba', 'Eu', 'Teff', + 'IMF1', 'IMF2', 'logfy', 'sigma2', 'velz2', 'logm7g', + 'hotteff', 'loghot', 'fy_logage', 'logtrans', 'logemline_h', + 'logemline_oiii', 'logemline_sii', 'logemline_ni', + 'logemline_nii', 'jitter', 'IMF3', 'logsky', 'IMF4', + 'h3', 'h4','ML_v','ML_i','ML_k','MW_v', 'MW_i','MW_k' + ]) + elif self.filetype == 50: + self.labels = np.array([ + 'chi2', 'velz', 'sigma', 'logage', 'zH', 'FeH', 'a', 'C', + 'N', 'Na', 'Mg', 'Si', 'K', 'Ca', 'Ti', 'V', 'Cr', 'Mn', + 'Co', 'Ni', 'Cu', 'Sr', 'Ba', 'Eu', 'Teff', 'IMF1', 'IMF2', + 'logfy', 'sigma2', 'velz2', 'logm7g', 'hotteff', 'loghot', + 'fy_logage', 'logtrans', 'logemline_h', 'logemline_oiii', + 'logemline_sii', 'logemline_ni', 'logemline_nii', 'jitter', + 'IMF3', 'logsky', 'IMF4', 'ML_v', 'ML_i', 'ML_k', 'MW_v', + 'MW_i', 'MW_k' + ]) results = Table(results, names=self.labels) @@ -72,16 +120,38 @@ def __init__(self, outfiles, read_mcmc=True, info=None, index=False): # is filled in abundance_correct() self.xFe = {} - self.results = results['Type', - 'chi2', 'velz', 'sigma', 'logage', 'zH', 'FeH', 'a', - 'C', 'N', 'Na', 'Mg', 'Si', 'K', 'Ca', 'Ti', 'V', 'Cr', - 'Mn', 'Co', 'Ni', 'Cu', 'Sr', 'Ba', 'Eu', 'Teff', - 'IMF1', 'IMF2', 'logfy', 'sigma2', 'velz2', 'logm7g', - 'hotteff', 'loghot', 'fy_logage', 'logemline_h', - 'logemline_oii', 'logemline_oiii', 'logemline_sii', - 'logemline_ni', 'logemline_nii', 'logtrans', 'jitter', - 'logsky', 'IMF3', 'IMF4', 'h3', 'h4', - 'ML_v','ML_i','ML_k','MW_v', 'MW_i','MW_k' ] + # Same filetype differentiations as above + if self.filetype == 53: + self.results = results['Type', + 'chi2', 'velz', 'sigma', 'logage', 'zH', 'FeH', 'a', + 'C', 'N', 'Na', 'Mg', 'Si', 'K', 'Ca', 'Ti', 'V', 'Cr', + 'Mn', 'Co', 'Ni', 'Cu', 'Sr', 'Ba', 'Eu', 'Teff', + 'IMF1', 'IMF2', 'logfy', 'sigma2', 'velz2', 'logm7g', + 'hotteff', 'loghot', 'fy_logage', 'logemline_h', + 'logemline_oii', 'logemline_oiii', 'logemline_sii', + 'logemline_ni', 'logemline_nii', 'logtrans', 'jitter', + 'logsky', 'IMF3', 'IMF4', 'h3', 'h4', + 'ML_v','ML_i','ML_k','MW_v', 'MW_i','MW_k' ] + elif self.filetype == 52: + self.results = results['Type', + 'chi2', 'velz', 'sigma', 'logage', 'zH', 'FeH', 'a', + 'C', 'N', 'Na', 'Mg', 'Si', 'K', 'Ca', 'Ti', 'V', 'Cr', + 'Mn', 'Co', 'Ni', 'Cu', 'Sr', 'Ba', 'Eu', 'Teff', + 'IMF1', 'IMF2', 'logfy', 'sigma2', 'velz2', 'logm7g', + 'hotteff', 'loghot', 'fy_logage', 'logtrans', 'logemline_h', + 'logemline_oiii', 'logemline_sii', 'logemline_ni', + 'logemline_nii', 'jitter', 'IMF3', 'logsky', 'IMF4', + 'h3', 'h4','ML_v','ML_i','ML_k','MW_v', 'MW_i','MW_k'] + elif self.filetype == 50: + self.results = results['Type', + 'chi2', 'velz', 'sigma', 'logage', 'zH', 'FeH', 'a', 'C', + 'N', 'Na', 'Mg', 'Si', 'K', 'Ca', 'Ti', 'V', 'Cr', 'Mn', + 'Co', 'Ni', 'Cu', 'Sr', 'Ba', 'Eu', 'Teff', 'IMF1', 'IMF2', + 'logfy', 'sigma2', 'velz2', 'logm7g', 'hotteff', 'loghot', + 'fy_logage', 'logtrans', 'logemline_h', 'logemline_oiii', + 'logemline_sii', 'logemline_ni', 'logemline_nii', 'jitter', + 'IMF3', 'logsky', 'IMF4', 'ML_v', 'ML_i', 'ML_k', 'MW_v', + 'MW_i', 'MW_k'] """ Read in input data and best fit model @@ -89,31 +159,28 @@ def __init__(self, outfiles, read_mcmc=True, info=None, index=False): This isn't going to work correctly if the file doesn't exist """ - #try: - m = np.loadtxt('{0}.bestspec'.format(self.outfiles)) - #except: - # warning = ('Do not have the *.bestspec file') - # warnings.warn(warning) - data = {} - data['wave'] = m[:,0]/(1.+self.results['velz'][5]*1e3/constants.c) - data['m_flux'] = m[:,1] # Model spectrum, normalization applied - data['d_flux'] = m[:,2] # Data spectrum - data['snr'] = m[:,3] # Including jitter and inflated errors - data['unc'] = 1/m[:,3] if not index: - data['poly'] = m[:,4] # Polynomial used to create m_flux - data['residual'] = (m[:,1] - m[:,2])/m[:,1] * 1e2 - self.spectra = data - - try: - m = np.loadtxt('{0}.bestspec2'.format(self.outfiles)) - model = {} - model['wave'] = m[:,0] - #model['wave'] = m[:,0]/(1.+self.results['velz'][5]*1e3/constants.c) - model['flux'] = m[:,1] - self.ext_model = model - except: - self.ext_model = None + m = np.loadtxt('{0}.bestspec'.format(self.outfiles)) + data = {} + data['wave'] = m[:,0]/(1.+self.results['velz'][5]*1e3/constants.c) + data['m_flux'] = m[:,1] # Model spectrum, normalization applied + data['d_flux'] = m[:,2] # Data spectrum + data['snr'] = m[:,3] # Including jitter and inflated errors + data['unc'] = 1/m[:,3] + if not index: + data['poly'] = m[:,4] # Polynomial used to create m_flux + data['residual'] = (m[:,1] - m[:,2])/m[:,1] * 1e2 + self.spectra = data + + try: + m = np.loadtxt('{0}.bestspec2'.format(self.outfiles)) + model = {} + model['wave'] = m[:,0] + #model['wave'] = m[:,0]/(1.+self.results['velz'][5]*1e3/constants.c) + model['flux'] = m[:,1] + self.ext_model = model + except: + self.ext_model = None """ Check the values of the nuisance parameters @@ -125,7 +192,7 @@ def __init__(self, outfiles, read_mcmc=True, info=None, index=False): # warnings.warn(warning.format(self.path, 'loghot', # self.results['loghot'][0])) - def get_total_met(self): + def get_total_met(self, return_posterior=False): zh = np.where(self.labels == 'zH') feh = np.where(self.labels == 'FeH') @@ -134,6 +201,9 @@ def get_total_met(self): #Computing errors directly from the chains. self.tmet = self.get_cls(total_met) + if return_posterior is True: + return total_met + def normalize_spectra(self): """ Normalize the data and model spectra @@ -146,6 +216,13 @@ def normalize_spectra(self): min_ = min(self.spectra['wave']) max_ = max(self.spectra['wave']) num = int((max_ - min_)/chunks) + 1 + + # If you have a model in .bestspec2, need to normalize this as well + if self.ext_model != None: + self.ext_model['flux'] = deepcopy(self.ext_model['flux']) + min_ext = min(self.ext_model['wave']) + max_ext = max(self.ext_model['wave']) + num_ext = int((max_ext - min_ext)/chunks) + 1 for i in range(num): k = ((self.spectra['wave'] >= min_ + chunks*i) & @@ -164,6 +241,18 @@ def normalize_spectra(self): self.spectra['m_flux_norm'][k], 2) poly = chebval(self.spectra['wave'][k], coeffs) self.spectra['m_flux_norm'][k] = self.spectra['m_flux_norm'][k]/poly + + # .bestspec2 model normalization + if self.ext_model != None: + for i in range(num_ext): + k_ext = ((self.ext_model['wave'] >= min_ext + chunks*i) & (self.ext_model['wave'] <= min_ext + chunks*(i+1))) + + if len(self.ext_model['flux'][k_ext]) < 10: + continue + + coeffs_ext = chebfit(self.ext_model['wave'][k_ext], self.ext_model['flux'][k_ext], 2) + poly_ext = chebval(self.ext_model['wave'][k_ext], coeffs_ext) + self.ext_model['flux'][k_ext] = self.ext_model['flux'][k_ext]/poly_ext def abundance_correct(self, s07=False, b14=False, m11=True): """ @@ -444,16 +533,138 @@ def plot_posterior(self, fname): def get_cls(self, distribution): distribution = np.sort(np.squeeze(distribution)) - num = self.nwalkers*self.nchain/self.nsample + if self.filetype == 50: + num = self.model_info['Nwalkers']*self.model_info['Nchain'] + else: + num = self.model_info['Nwalkers']*self.model_info['Nchain']/self.model_info['Nsample'] + ## ND: New file reading code, don't need this after rehaul. Delete. + #num = self.nwalkers*self.nchain/self.nsample lower = distribution[int(0.160*num)] median = distribution[int(0.500*num)] upper = distribution[int(0.840*num)] + cl95 = distribution[int(0.950*num)] std = np.std(distribution) - return {'cl50': median, 'cl84': upper, 'cl16': lower, 'std': std} + return {'cl50': median, 'cl84': upper, 'cl16': lower, + 'cl95': cl95, 'std': std} def write_params(self): pass + def get_mass(self): + """ + Based on the alf function. + + Compute amss in stars and remnants (normalized to 1 Msun at t=0) + for different IMF parameterizations (including a non-parameterization) + """ + + #try: # initial sanity check + # mlo < m2 + #except: + # print("Something is wrong: ml>m2") + + """ + mass = 0.0 + + if self.results['IMF4'] is None: + # normalize the weights so that 1 Msun formed at t=0 + imf_norm = (m2**(-self.results['IMF1']+2)-mlo**(-self.results['IMF1']+2))/(-self.results['IMF1']+2) + + m2**(-self.results['IMF1']+self.results['IMF2'])*(m3**(-self.results['IMF2']+2)- + m2**(-self.results['IMF2']+2))/(-self.['IMF2']+2) + + m2**(-self.['IMF1']+self.results['IMF2'])*(imf_hi**(-imf_up+2)-m3**(-imf_up+2))/(-imf_up+2) + + # mass from stars that are still alive + mass = (m2**(-self.results['IMF1']+2)-m_lo**(-self.results['IMF1']+2))/(-self.results['IMF1']+2) + if msto < m3: + mass = mass + m2**(-self.results['IMF1'+self.results['IMF2'])* + (msto**(-self.results['IMF2']+2)-m2**(-self.results['IMF2']+2))/(-self.['IMF2']+2) + else: + mass = mass + + m2**(-self.reslts['IMF1']+self.results['IMF2'])*(m3**(-self.results['IMF2']+2)-m2**(-self.results['IMF2']+2))/(-self.results['IMF2']+2) + + m2**(-self.results['IMF1']+self.results['IMF2'])*(msto**(-imf_up+2)-m3**(-imf_up+2))/(-imf_up+2) + + mass = mass/imf_norm + + + # Mass from blackhole remnants + # 40 < M < imf_up, leave behind a 0.5*M BH + mass = mass + 0.5*m2**(-self.results['IMF1']+self.results['IMF2'])*(imf_hi**(-imf_up+2)-bh_lim**(-imf_up+2))/(-imf_up+2)/imf_norm + + # Neutron star remnants + # 8.5 < M < 40, leave behind 1.4 MSun NS + mass = mass + 1.4*m2**(-self.results['IMF1']+self.results['IMF2'])*(bh_lim**(-imf_up+1)-ns_lim**(-imf_up+1))/(-imf_up+1)/imf_norm + + # White dwarf remnants + # M < 8.5, leave behind 0.077*M+0.58 WD + if msto < m3: + mass = mass + 0.48*m2**(-self.results['IMF1']+self.results['IMF2'])* + (ns_lim**(-imf_up+1)-m3**(-imf_up+1))/(-imf_up+1)/imf+norm + mass = mass + 0.48*m2**(-self.results['IMF1']+self.results['IMF2'])* + (m3**(-self.results['IMF2']+1)-msto**(-self.results['IMF2']+1))/(-self.results['IMF2']+1)/imf_norm + mass = mass + 0.077*m2**(-self.results['IMF1']+self.results['IMF2'])*(ns_lim**(-imf_up+2)-m3**(-imf_up+2))/(-imf_up+2)/imfnorm + mass = mass + + 0.077*m2**(-self.results['IMF1']+self.results['IMF2'])*(m3**(-self.results['IMF2']+2)-mtso**(-self.results['IMF2']+2))/(-self.results['IMF2']+2)/imf_norm + else: + mass = mass + 0.48*m2**(-self.results['IMF1']+self.results['IMF2'])*(ns_lim**(-imf_up+1)-mtso**(-imf_up+1))/(-imf_up+1)/imfnorm + mass = mass + 0.077*m2**(-self.results['IMF1']+self.results['IMF2'])*(ns_lim**(-imf_up+2)-msto**(-imf_up+2))/(-imf_up+2)/imfnorm + + + #else: Non-parametric IMF fit + + + #IF (PRESENT(timfnorm)) timfnorm = imfnorm + + return norm + """ + return None + + def get_m2l(self, filters): + """ + Based on the alf code. + + Computes mass-to-light ratios (AB mags) over user given filters + """ + + # convert to the "proper" units + aspec = self.data['d_flux']*lsun/1E6*self.data['wave']**2/clight/1E8/4/mypi/pc2cm**2 + + if self.model_info['mwimf'] == 1: # model was forced to fit for a MW (Kroupa) IMF + mass = get_mass(imf_lo, msto, kroupa_imf1, kroupa_imf2, kroupa_imf3) + """ + else: + if imf_type == 0: # single powerlaw IMF with a fixed lower-mass cutoff + mass = get_mass(imf_lo, msto, self.results['IMF1'], self.results['IMF1'], kroupa_imf3) + elif imf_type == 1: # double powerlaw with a fixed lower-mass cutoff + mass = get_mass(imf_lo, msto, self.results['IMF1'], self.results['IMF2'], kroupa_imf3) + elif imf_type == 2: # single powerlaw index with variable lower-mass cutoff + mass = get_mass(self.results['IMF3'], msto, self.results['IMF1'], self.results['IMF1'], kroupa_imf3) + elif imf_type == 3: # double powerlaw index with variable lower-mass cutoff + mass = get_mass(self.results['IMF3'], msto, self.results['IMF1'], self.results['IMF2'], kroupa_imf3) + elif imf_type == 4: # non-parametric IMF for 0.08-1.0 Msun; Salpeter slope at >1 Msun + mass = get_mass(imf_lo, msto, self.results['IMF1'], self.results['IMF2'], kroupa_imf3, + self.results['IMF3'], self.results['IMF4']) + """ + + + m2l = {} + for name, bandpass in filters: + mag = tsum(wave, aspec*f/lam) + if mag < 0.0: + m2l[name] == 0.0 + else: + mag = -2.5*np.log10(mag) - 48.60 + m2l[name] = mass//10**(2./5*(magsun[name]-mag)) + # if m2l > 100 then m2l = 0.0 + + return m2l + if __name__=='__main__': pass + + + + + + +