diff --git a/docs/parameters.rst b/docs/parameters.rst index 10b60ee..52f1070 100644 --- a/docs/parameters.rst +++ b/docs/parameters.rst @@ -136,6 +136,9 @@ An example config file can be generated by executing the following command (TBD) * - AMIN - str - (Optional) Minimum value for amplitude, in normal logarithmic scale. + * - DUST_MODEL + - int + - Dust attenuation model (default 0). 0 for Calzetti. 1 for MW. 2 for LCM. 3 for SMC. 4 for Kriek & Conroy (2013). * - AVMAX - float - (Optional) Maximum value for Av (dust attenuation in V-band), in mag. @@ -171,7 +174,7 @@ An example config file can be generated by executing the following command (TBD) - -**Parameters for a specific target:** +**Parameters specific for the target:** .. list-table:: :widths: 10 5 20 diff --git a/gsf/fitting.py b/gsf/fitting.py index bee457c..1d34656 100644 --- a/gsf/fitting.py +++ b/gsf/fitting.py @@ -16,12 +16,13 @@ import timeit from scipy import stats from scipy.stats import norm +import scipy.interpolate as interpolate from astropy.io import fits,ascii import corner # import from custom codes from .function import check_line_man, check_line_cz_man, calc_Dn4, savecpkl, get_leastsq -from .zfit import check_redshift +from .zfit import check_redshift,get_chi2 from .writing import get_param from .function_class import Func from .minimizer import Minimizer @@ -93,6 +94,7 @@ def update_input(self, inputs, c:float=3e18, Mpc_cm:float=3.08568025e+24, m0set: self.pixelscale = pixelscale self.Lsun = Lsun self.sigz = sigz + self.fitc_cz_prev = None # Magzp; try: @@ -183,12 +185,15 @@ def update_input(self, inputs, c:float=3e18, Mpc_cm:float=3.08568025e+24, m0set: # self.af = asdf.open(self.DIR_TMP + 'spec_all_' + self.ID + '_PA' + self.PA + '.asdf') try: self.DIR_EXTR = inputs['DIR_EXTR'] + except: + self.DIR_EXTR = False + + try: # Scaling for grism; self.Cz0 = float(inputs['CZ0']) self.Cz1 = float(inputs['CZ1']) self.Cz2 = float(inputs['CZ2']) except: - self.DIR_EXTR = False self.Cz0 = 1 self.Cz1 = 1 self.Cz2 = 1 @@ -227,7 +232,7 @@ def update_input(self, inputs, c:float=3e18, Mpc_cm:float=3.08568025e+24, m0set: self.logUMIN = self.logUFIX self.logUMAX = self.logUFIX self.DELlogU = 0 - self.logUs = np.arange(self.logUMIN, self.logUMAX, self.DELlogU) + self.logUs = np.asarray([self.logUMAX]) except: self.logUFIX = None else: @@ -405,6 +410,7 @@ def update_input(self, inputs, c:float=3e18, Mpc_cm:float=3.08568025e+24, m0set: self.has_ZFIX = True else: self.Zall = np.arange(self.Zmin, self.Zmax, self.delZ) + # If BPASS; if self.f_bpass == 1: try: @@ -417,6 +423,7 @@ def update_input(self, inputs, c:float=3e18, Mpc_cm:float=3.08568025e+24, m0set: self.Zsun = 0.020 Zbpass = [1e-5, 1e-4, 0.001, 0.002, 0.003, 0.004, 0.006, 0.008, 0.010, 0.020, 0.030, 0.040] Zbpass = np.log10(np.asarray(Zbpass)/self.Zsun) + try: # If ZFIX is found; iiz = np.argmin(np.abs(Zbpass[:] - float(inputs['ZFIX']) ) ) if Zbpass[iiz] - float(inputs['ZFIX']) != 0: @@ -437,7 +444,9 @@ def update_input(self, inputs, c:float=3e18, Mpc_cm:float=3.08568025e+24, m0set: self.delZ = 0.0001 self.Zmax,self.Zmin = np.max(self.Zall), np.min(self.Zall) print('Final list for log(Z_BPASS/Zsun) is:',self.Zall) - + if len(self.Zall)>1: + self.has_ZFIX = False + self.ZFIX = None # N of param: self.has_AVFIX = False @@ -799,7 +808,6 @@ def search_redshift(self, dict, xm_tmp, fm_tmp, zliml:float=0.01, zlimu:float=6. chi2s : numpy.array Array of chi2 values corresponding to zspace. ''' - import scipy.interpolate as interpolate zspace = np.arange(zliml,zlimu,delzz) chi2s = np.zeros((len(zspace),2), 'float') @@ -878,13 +886,19 @@ def residual_z(pars,z): def fit_redshift(self, xm_tmp, fm_tmp, delzz=0.01, ezmin=0.01, zliml=0.01, zlimu=None, snlim=0, priors=None, f_bb_zfit=True, f_line_check=False, - f_norm=True, f_lambda=False): + f_norm=True, f_lambda=False, zmax=20, fit_photometry=False, f_exclude_negative=False): ''' + Purpose + ------- Find the best-fit redshift, before going into a big fit, through an interactive inspection. This module is effective only when spec data is provided. + When spectrun is provided, this does redshift fit, but **by using the SED model guessed from the BB photometry**. + Thus, providing good BB photometry is critical in this step. Parameters ---------- + xm_tmp, fm_tmp : float array + SED model. delzz : float Delta z in redshift search space zliml : float @@ -901,35 +915,47 @@ def fit_redshift(self, xm_tmp, fm_tmp, delzz=0.01, ezmin=0.01, zliml=0.01, If True, line masking. priors : dict, optional Dictionary that contains z (redshift grid) and chi2 (chi-square). + f_exclude_negative : bool + Exclude negative fluxes in spectrum. Notes ----- Spectral data must be provided to make this work. ''' - import scipy.interpolate as interpolate + self.file_zprob = self.DIR_OUT + 'zprob_' + self.ID + '.txt' # NMC for zfit self.nmc_cz = int(self.inputs['NMCZ']) # For z prior. - #zliml = self.zgal - 0.5 - #if zlimu == None: - # zlimu = self.zgal + 0.5 - zliml = self.zmcmin - zlimu = self.zmcmax + if self.zmcmin != None: + zliml = self.zmcmin + else: + zliml = 0 + if self.zmcmax != None: + zlimu = self.zmcmax + else: + zlimu = 20 # Observed data; sn = self.dict['fy'] / self.dict['ey'] # Only spec data? - con_cz = (self.dict['NR']snlim) + if fit_photometry: + con_cz = ()#(sn>snlim) + else: + con_cz = (self.dict['NR']snlim) + if len(self.dict['fy'][con_cz])==0: if f_bb_zfit: - con_cz = (sn>snlim) + con_cz = ()#(sn>snlim) else: return 'y' + if f_exclude_negative: + con_cz &= (sn>snlim) + fy_cz = self.dict['fy'][con_cz] # Already scaled by self.Cz0 ey_cz = self.dict['ey'][con_cz] x_cz = self.dict['x'][con_cz] # Observed range @@ -953,12 +979,12 @@ def fit_redshift(self, xm_tmp, fm_tmp, delzz=0.01, ezmin=0.01, zliml=0.01, zprob = dprob[:,0] cprob = dprob[:,1] # Then interpolate to a common z grid; - zz_prob = np.arange(0,13,delzz) + zz_prob = np.arange(0,zmax,delzz) cprob_s = np.interp(zz_prob, zprob, cprob) prior_s = np.exp(-0.5 * cprob_s) prior_s /= np.sum(prior_s) else: - zz_prob = np.arange(0,13,delzz) + zz_prob = np.arange(zliml,zlimu,delzz) if priors != None: zprob = priors['z'] cprob = priors['chi2'] @@ -969,10 +995,7 @@ def fit_redshift(self, xm_tmp, fm_tmp, delzz=0.01, ezmin=0.01, zliml=0.01, prior_s[con_pri] = 0 if f_norm: prior_s /= np.sum(prior_s) - #prior_s /= np.sum(prior_s) - else: - zz_prob = np.arange(0,13,delzz) prior_s = zz_prob * 0 + 1. prior_s /= np.sum(prior_s) @@ -989,16 +1012,23 @@ def fit_redshift(self, xm_tmp, fm_tmp, delzz=0.01, ezmin=0.01, zliml=0.01, data_model[:,1] = fm_s data_model[:,2] = fy_cz data_model[:,3] = ey_cz - data_model_sort = np.sort(data_model, axis=0) + + data_model_sort = data_model[data_model[:, 0].argsort()] + #data_model_sort = np.sort(data_model, axis=0) # This does not work!! + plt.plot(data_model_sort[:,0], data_model_sort[:,1], 'gray', linestyle='--', linewidth=0.5, label='') # Model based on input z. plt.plot(data_model_sort[:,0], data_model_sort[:,2],'.b', linestyle='-', linewidth=0.5, label='Obs.') # Observation plt.errorbar(data_model_sort[:,0], data_model_sort[:,2], yerr=data_model_sort[:,3], color='b', capsize=0, linewidth=0.5) # Observation + # Write prob distribution; + if True: + get_chi2(zz_prob, fy_cz, ey_cz, x_cz, fm_tmp, xm_tmp/(1+self.zgal), self.file_zprob) + print('############################') print('Start MCMC for redshift fit') print('############################') res_cz, fitc_cz = check_redshift(fy_cz, ey_cz, x_cz, fm_tmp, xm_tmp/(1+self.zgal), self.zgal, self.z_prior, self.p_prior, \ - NR_cz, zliml, zlimu, self.nmc_cz, self.nwalk_cz) + NR_cz, zliml, zlimu, self.nmc_cz, self.nwalk_cz, include_photometry=False) z_cz = np.percentile(res_cz.flatchain['z'], [16,50,84]) scl_cz0 = np.percentile(res_cz.flatchain['Cz0'], [16,50,84]) scl_cz1 = np.percentile(res_cz.flatchain['Cz1'], [16,50,84]) @@ -1019,6 +1049,7 @@ def fit_redshift(self, xm_tmp, fm_tmp, delzz=0.01, ezmin=0.01, zliml=0.01, print('Recommended redshift, Cz0, Cz1, and Cz2, %.5f %.5f %.5f %.5f, with chi2/nu=%.3f'%(zrecom, Czrec0, Czrec1, Czrec2, fitc_cz[1])) print('\n\n') fit_label = 'Proposed model' + #self.fitc_cz = fitc_cz[1] else: print('fzvis is set to False. z fit not happening.') @@ -1032,7 +1063,7 @@ def fit_redshift(self, xm_tmp, fm_tmp, delzz=0.01, ezmin=0.01, zliml=0.01, print('Redshift error is assumed to %.1f.'%(ezmin)) z_cz = [zmcmin, self.zprev, zmcmax] - zrecom = z_cz[1] + zrecom = z_cz[1] scl_cz0 = [1.,1.,1.] scl_cz1 = [1.,1.,1.] scl_cz2 = [1.,1.,1.] @@ -1040,6 +1071,7 @@ def fit_redshift(self, xm_tmp, fm_tmp, delzz=0.01, ezmin=0.01, zliml=0.01, Czrec1 = scl_cz1[1] Czrec2 = scl_cz2[1] res_cz = None + fitc_cz = [99,99] # If this label is being used, it means that the fit is failed. fit_label = 'Current model' @@ -1067,7 +1099,8 @@ def fit_redshift(self, xm_tmp, fm_tmp, delzz=0.01, ezmin=0.01, zliml=0.01, data_model_new = np.zeros((len(x_cz),4),'float') data_model_new[:,0] = x_cz data_model_new[:,1] = fm_s - data_model_new_sort = np.sort(data_model_new, axis=0) + # data_model_new_sort = np.sort(data_model_new, axis=0) + data_model_new_sort = data_model_new[data_model_new[:, 0].argsort()] plt.plot(data_model_new_sort[:,0], data_model_new_sort[:,1], 'r', linestyle='-', linewidth=0.5, label='%s ($z=%.5f$)'%(fit_label,zrecom)) # Model based on recomended z. plt.plot(x_cz[con_line], fm_s[con_line], color='orange', marker='o', linestyle='', linewidth=3.) @@ -1092,14 +1125,16 @@ def fit_redshift(self, xm_tmp, fm_tmp, delzz=0.01, ezmin=0.01, zliml=0.01, data_obsbb[:,0],data_obsbb[:,1] = self.dict['xbb'],self.dict['fybb'] if len(fm_tmp) == len(self.dict['xbb']): # BB only; data_obsbb[:,2] = fm_tmp - data_obsbb_sort = np.sort(data_obsbb, axis=0) - + #data_obsbb_sort = np.sort(data_obsbb, axis=0) + data_obsbb_sort = data_obsbb[data_obsbb[:, 0].argsort()] + if len(fm_tmp) == len(self.dict['xbb']): # BB only; plt.scatter(data_obsbb_sort[:,0], data_obsbb_sort[:,2], color='none', marker='d', s=50, edgecolor='gray', zorder=4, label='Current model ($z=%.5f$)'%(self.zgal)) else: model_spec = np.zeros((len(fm_tmp),2), 'float') model_spec[:,0],model_spec[:,1] = xm_tmp,fm_tmp - model_spec_sort = np.sort(model_spec, axis=0) + #model_spec_sort = np.sort(model_spec, axis=0) + model_spec_sort = model_spec[model_spec[:, 0].argsort()] plt.plot(model_spec_sort[:,0], model_spec_sort[:,1], marker='.', color='gray', ms=1, linestyle='-', linewidth=0.5, zorder=4, label='Current model ($z=%.5f$)'%(self.zgal)) try: @@ -1138,11 +1173,16 @@ def fit_redshift(self, xm_tmp, fm_tmp, delzz=0.01, ezmin=0.01, zliml=0.01, print('##############################################################\n') plt.show() - flag_z = raw_input('Do you want to continue with the input redshift, Cz0, Cz1 and Cz2, %.5f %.5f %.5f %.5f? ([y]/n/m) '%\ - (self.zgal, self.Cz0, self.Cz1, self.Cz2)) + flag_z = raw_input('Do you want to continue with the input redshift, Cz0, Cz1, Cz2, and chi2/nu, %.5f %.5f %.5f %.5f %.5f? ([y]/n/m) '%\ + (self.zgal, self.Cz0, self.Cz1, self.Cz2, self.fitc_cz_prev)) else: flag_z = 'y' + try: + self.z_cz_prev = self.z_cz + except: + self.z_cz_prev = [self.zgal,self.zgal,self.zgal] + # Write it to self; self.zrecom = zrecom self.Czrec0 = Czrec0 * self.Cz0 @@ -1153,12 +1193,15 @@ def fit_redshift(self, xm_tmp, fm_tmp, delzz=0.01, ezmin=0.01, zliml=0.01, self.scl_cz1 = scl_cz1 self.scl_cz2 = scl_cz2 self.res_cz = res_cz + self.fitc_cz_prev = fitc_cz[1] return flag_z def get_zdist(self, f_interact=False, f_ascii=True): ''' + Purpose + ------- Saves a plot of z-distribution. Parameters @@ -1180,11 +1223,11 @@ def get_zdist(self, f_interact=False, f_ascii=True): (self.z_cz[1],self.z_cz[1]-self.z_cz[0],self.z_cz[2]-self.z_cz[1], self.Cz0, self.Cz1, self.Cz2)) if f_ascii: - file_ascii_out = self.DIR_OUT + 'zprob_' + self.ID + '.txt' + file_ascii_out = self.DIR_OUT + 'zmc_' + self.ID + '.txt' fw_ascii = open(file_ascii_out,'w') fw_ascii.write('# z pz\n') - for ii in range(len(xx)): - fw_ascii.write('%.3f %.3f\n'%(xx[ii],yy[ii])) + for ii in range(len(n)): + fw_ascii.write('%.3f %.3f\n'%(nbins[ii]+(nbins[ii+1]-nbins[ii])/2.,n[ii])) fw_ascii.close() xx = yy * 0 + self.z_cz[0] @@ -1205,21 +1248,46 @@ def get_zdist(self, f_interact=False, f_ascii=True): ax1.legend(loc=0) # Save: - file_out = self.DIR_OUT + 'zprob_' + self.ID + '.png' + file_out = self.DIR_OUT + 'zmc_' + self.ID + '.png' print('Figure is saved in %s'%file_out) if f_interact: fig.savefig(file_out, dpi=300) - return fig, ax1 + # return fig, ax1 else: plt.savefig(file_out, dpi=300) plt.close() - return True + # return True except: print('z-distribution figure is not generated.') pass + if os.path.exists(self.file_zprob): + # Also, make zprob; + fig = plt.figure(figsize=(6.5,2.5)) + fig.subplots_adjust(top=0.96, bottom=0.16, left=0.09, right=0.99, hspace=0.15, wspace=0.25) + ax1 = fig.add_subplot(111) + + fd = ascii.read(self.file_zprob) + z = fd['z'] + pz = fd['p(z)'] + pz /= pz.max() + + # prob: + ax1.plot(z, pz, linestyle='-', linewidth=1, color='r', label='') + + # Label: + ax1.set_xlabel('Redshift') + ax1.set_ylabel('$p(z)$') + file_out = self.DIR_OUT + 'zprob_' + self.ID + '.png' + if f_interact: + fig.savefig(file_out, dpi=300) + else: + plt.savefig(file_out, dpi=300) + plt.close() + + def add_param(self, fit_params, sigz=1.0, zmin=None, zmax=None): ''' Add new parameters. @@ -1228,10 +1296,10 @@ def add_param(self, fit_params, sigz=1.0, zmin=None, zmax=None): f_add = False # Redshift if self.fzmc == 1: - if zmin == None: - self.zmcmin = self.zgal-(self.z_cz[1]-self.z_cz[0])*sigz - if zmax == None: - self.zmcmax = self.zgal+(self.z_cz[2]-self.z_cz[1])*sigz + # if zmin == None: + # self.zmcmin = self.zgal-(self.z_cz[1]-self.z_cz[0])*sigz + # if zmax == None: + # self.zmcmax = self.zgal+(self.z_cz[2]-self.z_cz[1])*sigz fit_params.add('zmc', value=self.zgal, min=self.zmcmin, max=self.zmcmax) print('Redshift is set as a free parameter (z in [%.2f:%.2f])'%(self.zmcmin, self.zmcmax)) f_add = True @@ -1389,10 +1457,10 @@ def set_param(self): else: fit_params.add('Z'+str(aa), value=0, min=self.Zmin, max=self.Zmax) else: - try: + if self.has_ZFIX: aa = 0 fit_params.add('Z'+str(aa), value=self.ZFIX, vary=False) - except: + else: aa = 0 if np.min(self.Zall)==np.max(self.Zall): fit_params.add('Z'+str(aa), value=np.min(self.Zall), vary=False) @@ -1546,7 +1614,8 @@ def get_shuffle(self, out, nshuf=3.0, amp=1e-4): def main(self, cornerplot:bool=True, specplot=1, sigz=1.0, ezmin=0.01, ferr=0, f_move:bool=False, verbose:bool=False, skip_fitz:bool=False, out=None, f_plot_accept:bool=True, - f_shuffle:bool=True, amp_shuffle=1e-2, check_converge:bool=True, Zini=None, f_plot_chain:bool=True): + f_shuffle:bool=True, amp_shuffle=1e-2, check_converge:bool=True, Zini=None, f_plot_chain:bool=True, + f_chind:bool=True, ncpu:int=0): ''' Main module of this script. @@ -1600,15 +1669,19 @@ def main(self, cornerplot:bool=True, specplot=1, sigz=1.0, ezmin=0.01, ferr=0, #################################### # Get initial parameters if not skip_fitz or out == None: - + + # Do you want to prepare a template for redshift fit by using only spectrum?; + f_only_spec = False out, chidef, Zbest = get_leastsq(self, Zini, self.fneld, self.age, self.fit_params, class_post.residual,\ - self.dict['fy'], self.dict['ey'], self.dict['wht2'], self.ID) + self.dict['fy'], self.dict['ey'], self.dict['wht2'], self.ID, f_only_spec=f_only_spec) # Best fit csq = out.chisqr rcsq = out.redchi fitc = [csq, rcsq] # Chi2, Reduced-chi2 ZZ = Zbest # This is really important/does affect lnprob/residual. + if self.fitc_cz_prev == None: + self.fitc_cz_prev = rcsq print('\n\n') print('#####################################') @@ -1632,16 +1705,13 @@ def main(self, cornerplot:bool=True, specplot=1, sigz=1.0, ezmin=0.01, ferr=0, if skip_fitz: flag_z = 'y' else: - #con_bb = (nrd_tmp>=1e4) - flag_z = self.fit_redshift(xm_tmp, fm_tmp) + flag_z = self.fit_redshift(xm_tmp, fm_tmp, delzz=0.001, fit_photometry=False) ################################################# # Gor for mcmc phase ################################################# if flag_z == 'y' or flag_z == '': - self.get_zdist() - ####################### # Add parameters; ####################### @@ -1660,8 +1730,9 @@ def main(self, cornerplot:bool=True, specplot=1, sigz=1.0, ezmin=0.01, ferr=0, fit_name = self.fneld out = minimize(class_post.residual, self.fit_params, args=(self.dict['fy'], self.dict['ey'], self.dict['wht2'], self.f_dust), method=fit_name) - print('\nMinimizer refinement;') - print(fit_report(out)) + # showing this is confusing. + # print('\nMinimizer refinement;') + # print(fit_report(out)) # Fix params to what we had before. if self.fzmc: @@ -1685,7 +1756,6 @@ def main(self, cornerplot:bool=True, specplot=1, sigz=1.0, ezmin=0.01, ferr=0, ################################ print('\nMinimizer Defined\n') - ncpu = 0 print('########################') print('### Starting sampling ##') @@ -1716,7 +1786,7 @@ def main(self, cornerplot:bool=True, specplot=1, sigz=1.0, ezmin=0.01, ferr=0, sampler = zeus.EnsembleSampler(self.nwalk, self.ndim, class_post.lnprob_emcee, \ args=[out.params, self.dict['fy'], self.dict['ey'], self.dict['wht2'], self.f_dust], \ moves=moves, maxiter=1e6,\ - kwargs={'f_val':True, 'out':out, 'lnpreject':-np.inf},\ + kwargs={'f_val':True, 'out':out, 'lnpreject':-np.inf, 'f_chind':f_chind},\ ) # Run MCMC nburn = int(self.nmc/10) @@ -2004,8 +2074,12 @@ def __init__(self, flatchain, var_names, params_value, res): else: print('\n\n') - flag_gen = raw_input('Do you want to make templates with recommended redshift, Cz0, Cz1, and Cz2 , %.5f %.5f %.5f %.5f? ([y]/n) '%\ - (self.zrecom, self.Czrec0, self.Czrec1, self.Czrec2)) + + flag_gen = raw_input('Do you want to make templates with recommended redshift, Cz0, Cz1, Cz2, and chi2/nu , %.5f %.5f %.5f %.5f %.5f? ([y]/n) '%\ + (self.zrecom, self.Czrec0, self.Czrec1, self.Czrec2, self.fitc_cz_prev)) + + self.get_zdist() + if flag_gen == 'y' or flag_gen == '': self.zprev = self.zgal # Input redshift for previous run self.zgal = self.zrecom # Recommended redshift from previous run diff --git a/gsf/function.py b/gsf/function.py index d02a85f..83212cd 100644 --- a/gsf/function.py +++ b/gsf/function.py @@ -248,7 +248,8 @@ def loadcpkl(cpklfile): return data -def get_leastsq(MB, ZZtmp, fneld, age, fit_params, residual, fy, ey, wht, ID0, chidef=None, Zbest=0, f_keep=False): +def get_leastsq(MB, ZZtmp, fneld, age, fit_params, residual, fy, ey, wht, ID0, + chidef=None, Zbest=0, f_keep=False, f_only_spec=False): ''' Get initial parameters at various Z ''' @@ -283,7 +284,10 @@ def get_leastsq(MB, ZZtmp, fneld, age, fit_params, residual, fy, ey, wht, ID0, c fit_params['Z'+str(aa)].value = ZZ f_fir = False - out_tmp = minimize(residual, fit_params, args=(fy, ey, wht, f_fir), method=fit_name) # nelder is the most efficient. + + out_tmp = minimize(residual, fit_params, args=(fy, ey, wht, f_fir), + method=fit_name, kws={'f_only_spec':f_only_spec}) + csq = out_tmp.chisqr rcsq = out_tmp.redchi fitc = [csq, rcsq] # Chi2, Reduced-chi2 @@ -550,12 +554,18 @@ def flamtonu(lam, flam, m0set=25.0): return fnu -def fnutolam(lam, fnu, m0set=25.0): +def fnutolam(lam, fnu, m0set=25.0, m0=-48.6): ''' Converts from Fnu to Flam, from mag zeropoint of m0set (to -48.6). - + + Parameters + ---------- + m0set : float + current magzp. + m0 : float + target magzp. The default, -48.6, is for flam (erg/s/cm2/lambda). ''' - Ctmp = lam**2/c * 10**((48.6+m0set)/2.5) + Ctmp = lam**2/c * 10**((m0set-m0)/2.5) flam = fnu / Ctmp return flam diff --git a/gsf/function_class.py b/gsf/function_class.py index 916953e..ec11b4d 100644 --- a/gsf/function_class.py +++ b/gsf/function_class.py @@ -224,7 +224,7 @@ def open_spec_fits_dir(self, nage:int, nz:int, kk, Av00:float, zgal:float, A00:f elif self.dust_model == 4: # Kriek&Conroy with gamma=-0.2 yyd, xxd, nrd = dust_kc(xx, yy, Av00, nr, Rv=4.05, gamma=-0.2) else: - print('No entry. Dust model is set to Calzetti') + print('No valid entry. Dust model is set to Calzetti') yyd, xxd, nrd = dust_calz(xx, yy, Av00, nr) xxd *= (1.+zgal) @@ -940,7 +940,7 @@ def open_spec_fits_dir(self, nage, nz, kk, Av00, zgal, A00): elif self.dust_model == 4: # Kriek&Conroy with gamma=-0.2 yyd, xxd, nrd = dust_kc(xx, yy, Av00, nr, Rv=4.05, gamma=-0.2) else: - print('No entry. Dust model is set to Calzetti') + print('No valid entry. Dust model is set to Calzetti') yyd, xxd, nrd = dust_calz(xx, yy, Av00, nr) xxd *= (1.+zgal) diff --git a/gsf/maketmp_filt.py b/gsf/maketmp_filt.py index a24960f..d23ab07 100644 --- a/gsf/maketmp_filt.py +++ b/gsf/maketmp_filt.py @@ -435,16 +435,23 @@ def maketemp(MB, ebblim=1e10, lamliml=0., lamlimu=50000., ncolbb=10000, # READ BB photometry from CAT_BB: ############################# if CAT_BB: + key_id = 'id' fd0 = ascii.read(CAT_BB) - id0 = fd0['id'].astype('str') + try: + id0 = fd0[key_id].astype('str') + except: + key_id = 'ID' + id0 = fd0[key_id].astype('str') + ii0 = np.where(id0[:]==ID) try: if len(ii0[0]) == 0: print('Could not find the column for [ID: %s] in the input BB catalog! Exiting.'%(ID)) return False - id = fd0['id'][ii0] + id = fd0[key_id][ii0] except: print('Could not find the column for [ID: %s] in the input BB catalog! Exiting.'%(ID)) + print(fd0) return False fbb = np.zeros(len(SFILT), dtype='float') @@ -494,18 +501,25 @@ def maketemp(MB, ebblim=1e10, lamliml=0., lamlimu=50000., ncolbb=10000, # Dust flux; if MB.f_dust: + key_id = 'id' fdd = ascii.read(CAT_BB_DUST) try: - id0 = fdd['id'].astype('str') + try: + id0 = fdd[key_id].astype('str') + except: + key_id = 'ID' + id0 = fdd[key_id].astype('str') + ii0 = np.where(id0[:]==ID) try: - id = fd0['id'][ii0] + id = fd0[key_id][ii0] except: print('Could not find the column for [ID: %s] in the input BB catalog! Exiting.'%(ID)) return False except: return False - id = fdd['id'] + + id = fdd[key_id] fbb_d = np.zeros(len(DFILT), dtype='float') ebb_d = np.zeros(len(DFILT), dtype='float') @@ -617,7 +631,7 @@ def maketemp(MB, ebblim=1e10, lamliml=0., lamlimu=50000., ncolbb=10000, try: spec_mul_nu_conv[ss,:] = convolve(spec_mul_nu[ss], LSF, boundary='extend') except: - print('Error. No convolution is happening...') + #print('Error. No convolution is happening...') spec_mul_nu_conv[ss,:] = spec_mul_nu[ss] if zz==0 and ss==0: print('Kernel is too small. No convolution.') diff --git a/gsf/maketmp_z0.py b/gsf/maketmp_z0.py index 3dd9523..5ae0682 100644 --- a/gsf/maketmp_z0.py +++ b/gsf/maketmp_z0.py @@ -9,7 +9,7 @@ INDICES = ['G4300', 'Mgb', 'Fe5270', 'Fe5335', 'NaD', 'Hb', 'Fe4668', 'Fe5015', 'Fe5709', 'Fe5782', 'Mg1', 'Mg2', 'TiO1', 'TiO2'] -def make_tmp_z0(MB, lammin=100, lammax=160000, tau_lim=0.001, force_no_neb=False, Zforce=None): +def make_tmp_z0(MB, lammin=100, lammax=160000, tau_lim=0.001, force_no_neb=False, Zforce=None, f_mp=True): ''' This is for the preparation of default template, with FSPS, at z=0. Should be run before SED fitting. diff --git a/gsf/plot_sed.py b/gsf/plot_sed.py index 8d28f98..d0e867a 100644 --- a/gsf/plot_sed.py +++ b/gsf/plot_sed.py @@ -17,6 +17,7 @@ col = ['violet', 'indigo', 'b', 'lightblue', 'lightgreen', 'g', 'orange', 'coral', 'r', 'darkred']#, 'k'] + def plot_sed(MB, flim=0.01, fil_path='./', scale=1e-19, f_chind=True, figpdf=False, save_sed=True, mmax=300, dust_model=0, DIR_TMP='./templates/', f_label=False, f_bbbox=False, verbose=False, f_silence=True, f_fill=False, f_fancyplot=False, f_Alog=True, dpi=300, f_plot_filter=True, f_plot_resid=False, NRbb_lim=10000): @@ -98,7 +99,6 @@ def gaus(x,a,x0,sigma): chimax = 1. m0set = MB.m0set Mpc_cm = MB.Mpc_cm - d = MB.d * scale ################## # Fitting Results @@ -289,14 +289,20 @@ def gaus(x,a,x0,sigma): AAAA BBBB """ - fig,axes = plt.subplot_mosaic(mosaic=fig_mosaic, figsize=(5.5,4.2)) - fig.subplots_adjust(top=0.98, bottom=0.16, left=0.1, right=0.99, hspace=0.15, wspace=0.25) + fig,axes = plt.subplot_mosaic(mosaic=fig_mosaic, figsize=(5.5,4.)) + fig.subplots_adjust(top=0.98, bottom=0.16, left=0.08, right=0.99, hspace=0.15, wspace=0.25) ax1 = axes['A'] else: - fig = plt.figure(figsize=(5.5,2.2)) - fig.subplots_adjust(top=0.98, bottom=0.16, left=0.1, right=0.99, hspace=0.15, wspace=0.25) + fig = plt.figure(figsize=(5.5,2.)) + fig.subplots_adjust(top=0.98, bottom=0.16, left=0.08, right=0.99, hspace=0.15, wspace=0.25) ax1 = fig.add_subplot(111) + # Determine scale here; + if scale == None: + conbb_hs = (fybb/eybb > SNlim) + scale = 10**(int(np.log10(np.nanmax(fybb[conbb_hs] * c / np.square(xbb[conbb_hs])) / MB.d))) / 10 + d = MB.d * scale + ####################################### # D.Kelson like Box for BB photometry ####################################### @@ -495,8 +501,8 @@ def gaus(x,a,x0,sigma): xboxl = 17000 xboxu = 28000 - ax1.set_xlabel('Observed wavelength ($\mathrm{\mu m}$)', fontsize=12) - ax1.set_ylabel('Flux ($10^{%d}\mathrm{erg}/\mathrm{s}/\mathrm{cm}^{2}/\mathrm{\AA}$)'%(np.log10(scale)),fontsize=12,labelpad=-2) + ax1.set_xlabel('Observed wavelength ($\mathrm{\mu m}$)', fontsize=11) + ax1.set_ylabel('Flux ($10^{%d}\mathrm{erg}/\mathrm{s}/\mathrm{cm}^{2}/\mathrm{\AA}$)'%(np.log10(scale)),fontsize=11,labelpad=-2) x1min = 2000 x1max = 100000 @@ -778,17 +784,15 @@ def gaus(x,a,x0,sigma): from astropy.convolution import convolve from .maketmp_filt import get_LSF LSF, _ = get_LSF(MB.inputs, MB.DIR_EXTR, ID, x1_tot[:]/(1.+zbes), c=3e18) - spec_grsm16 = convolve(ytmp16[:], LSF, boundary='extend') - spec_grsm50 = convolve(ytmp50[:], LSF, boundary='extend') - spec_grsm84 = convolve(ytmp84[:], LSF, boundary='extend') - ''' - plt.close() - plt.plot(x1_tot, ytmp50, color='r') - plt.plot(x1_tot, spec_grsm50, color='gray') - print(spec_grsm50) - plt.xlim(8000,20000) - plt.show() - ''' + try: + spec_grsm16 = convolve(ytmp16[:], LSF, boundary='extend') + spec_grsm50 = convolve(ytmp50[:], LSF, boundary='extend') + spec_grsm84 = convolve(ytmp84[:], LSF, boundary='extend') + except: + spec_grsm16 = ytmp16[:] + spec_grsm50 = ytmp50[:] + spec_grsm84 = ytmp84[:] + if True: ax2t.plot(x1_tot[:], ytmp50, '-', lw=0.5, color='gray', zorder=3., alpha=1.0) else: @@ -1012,11 +1016,11 @@ def func_tmp(xint,eobs,fmodel): # Grism; if f_grsm: - col2 = fits.Column(name='f_model_conv_16', format='E', unit='1e%derg/s/cm2/AA'%(np.log10(scale)), array=spec_grsm16) + col2 = fits.Column(name='f_model_conv_16', format='E', unit='1e%derg/s/cm2/AA'%(np.log10(scale)), array=spec_grsm16) col00.append(col2) - col3 = fits.Column(name='f_model_conv_50', format='E', unit='1e%derg/s/cm2/AA'%(np.log10(scale)), array=spec_grsm50) + col3 = fits.Column(name='f_model_conv_50', format='E', unit='1e%derg/s/cm2/AA'%(np.log10(scale)), array=spec_grsm50) col00.append(col3) - col4 = fits.Column(name='f_model_conv_84', format='E', unit='1e%derg/s/cm2/AA'%(np.log10(scale)), array=spec_grsm84) + col4 = fits.Column(name='f_model_conv_84', format='E', unit='1e%derg/s/cm2/AA'%(np.log10(scale)), array=spec_grsm84) col00.append(col4) # BB for dust @@ -1055,10 +1059,10 @@ def func_tmp(xint,eobs,fmodel): try: # Muv - MUV = -2.5 * np.log10(Fuv[:]) + 25.0 - hdr['MUV16'] = -2.5 * np.log10(np.percentile(Fuv[:],16)) + 25.0 - hdr['MUV50'] = -2.5 * np.log10(np.percentile(Fuv[:],50)) + 25.0 - hdr['MUV84'] = -2.5 * np.log10(np.percentile(Fuv[:],84)) + 25.0 + MUV = -2.5 * np.log10(Fuv[:]) + MB.m0set + hdr['MUV16'] = -2.5 * np.log10(np.percentile(Fuv[:],16)) + MB.m0set + hdr['MUV50'] = -2.5 * np.log10(np.percentile(Fuv[:],50)) + MB.m0set + hdr['MUV84'] = -2.5 * np.log10(np.percentile(Fuv[:],84)) + MB.m0set # Fuv (!= flux of Muv) hdr['FUV16'] = np.percentile(Fuv28[:],16) @@ -1198,12 +1202,12 @@ def func_tmp(xint,eobs,fmodel): if f_label: fd = fits.open(MB.DIR_OUT + 'SFH_' + ID + '.fits')[0].header if MB.f_dust: - label = 'ID: %s\n$z:%.2f$\n$\log M_\mathrm{*}/M_\odot:%.2f$\n$\log M_\mathrm{dust}/M_\odot:%.2f$\n$T_\mathrm{dust}/K:%.1f$\n$\log Z_\mathrm{*}/Z_\odot:%.2f$\n$\log T_\mathrm{*}$/Gyr$:%.2f$\n$A_V$/mag$:%.2f$\n$\\chi^2/\\nu:%.2f$'\ - %(ID, zbes, float(fd['Mstel_50']), MD50, TD50, float(fd['Z_MW_50']), float(fd['T_MW_50']), float(fd['AV_50']), fin_chi2) + label = 'ID: %s\n$z:%.2f$\n$\log M_\mathrm{*}/M_\odot:%.2f$\n$\log M_\mathrm{dust}/M_\odot:%.2f$\n$T_\mathrm{dust}/K:%.1f$\n$\log Z_\mathrm{*}/Z_\odot:%.2f$\n$\log T_\mathrm{*}$/Gyr$:%.2f$\n$A_V$/mag$:%.2f$'\ + %(ID, zbes, float(fd['Mstel_50']), MD50, TD50, float(fd['Z_MW_50']), float(fd['T_MW_50']), float(fd['AV_50']))#, fin_chi2) ylabel = ymax*0.45 else: - label = 'ID: %s\n$z:%.2f$\n$\log M_\mathrm{*}/M_\odot:%.2f$\n$\log Z_\mathrm{*}/Z_\odot:%.2f$\n$\log T_\mathrm{*}$/Gyr$:%.2f$\n$A_V$/mag$:%.2f$\n$\\chi^2/\\nu:%.2f$'\ - %(ID, zbes, float(fd['Mstel_50']), float(fd['Z_MW_50']), float(fd['T_MW_50']), float(fd['AV_50']), fin_chi2) + label = 'ID: %s\n$z:%.2f$\n$\log M_\mathrm{*}/M_\odot:%.2f$\n$\log Z_\mathrm{*}/Z_\odot:%.2f$\n$\log T_\mathrm{*}$/Gyr$:%.2f$\n$A_V$/mag$:%.2f$'\ + %(ID, zbes, float(fd['Mstel_50']), float(fd['Z_MW_50']), float(fd['T_MW_50']), float(fd['AV_50'])) ylabel = ymax*0.25 if f_grsm: @@ -1402,7 +1406,6 @@ def gaus(x,a,x0,sigma): chimax = 1. m0set = MB.m0set Mpc_cm = MB.Mpc_cm - d = MB.d * scale ################## # Fitting Results @@ -1614,6 +1617,12 @@ def gaus(x,a,x0,sigma): fig.subplots_adjust(top=0.98, bottom=0.16, left=0.1, right=0.99, hspace=0.15, wspace=0.25) ax1 = fig.add_subplot(111) + # Determine scale here; + if scale == None: + conbb_hs = (fybb/eybb > SNlim) + scale = 10**(int(np.log10(np.nanmax(fybb[conbb_hs] * c / np.square(xbb[conbb_hs])) / MB.d))) + d = MB.d * scale + ####################################### # D.Kelson like Box for BB photometry ####################################### @@ -2311,10 +2320,10 @@ def func_tmp(xint,eobs,fmodel): try: # Muv - MUV = -2.5 * np.log10(Fuv[:]) + 25.0 - hdr['MUV16'] = -2.5 * np.log10(np.percentile(Fuv[:],16)) + 25.0 - hdr['MUV50'] = -2.5 * np.log10(np.percentile(Fuv[:],50)) + 25.0 - hdr['MUV84'] = -2.5 * np.log10(np.percentile(Fuv[:],84)) + 25.0 + MUV = -2.5 * np.log10(Fuv[:]) + MB.m0set + hdr['MUV16'] = -2.5 * np.log10(np.percentile(Fuv[:],16)) + MB.m0set + hdr['MUV50'] = -2.5 * np.log10(np.percentile(Fuv[:],50)) + MB.m0set + hdr['MUV84'] = -2.5 * np.log10(np.percentile(Fuv[:],84)) + MB.m0set # Fuv (!= flux of Muv) hdr['FUV16'] = np.percentile(Fuv28[:],16) @@ -2621,7 +2630,6 @@ def plot_corner_physparam_summary(MB, fig=None, out_ind=0, DIR_OUT='./', mmax=10 ID = MB.ID Z = MB.Zall age = MB.age - d = MB.d * scale c = MB.c tau0 = MB.tau0 @@ -2762,6 +2770,12 @@ def plot_corner_physparam_summary(MB, fig=None, out_ind=0, DIR_OUT='./', mmax=10 ey = np.append(ey01,eg_bb) wht = 1./np.square(ey) + # Determine scale here; + if scale == None: + conbb_hs = (fybb/eybb > SNlim) + scale = 10**(int(np.log10(np.nanmax(fybb[conbb_hs] * c / np.square(xbb[conbb_hs])) / MB.d))) + d = MB.d * scale + # BB photometry conspec = (NR0) con_res_r = (wht==0) diff --git a/gsf/writing.py b/gsf/writing.py index 5c78752..26be777 100644 --- a/gsf/writing.py +++ b/gsf/writing.py @@ -24,7 +24,7 @@ def get_param(self, res, fitc, tcalc=1., burnin=-1): Czrec2 = self.Cz2 try: - z_cz = self.z_cz + z_cz = self.z_cz_prev scl_cz0 = self.scl_cz0 scl_cz1 = self.scl_cz1 scl_cz2 = self.scl_cz2 diff --git a/gsf/zfit.py b/gsf/zfit.py index 7086d32..64271b8 100644 --- a/gsf/zfit.py +++ b/gsf/zfit.py @@ -1,6 +1,96 @@ +import numpy as np +import sys +from lmfit import Model, Parameters, minimize, fit_report, Minimizer +import scipy.interpolate as interpolate +import matplotlib.pyplot as plt +from .function import check_line_cz_man + + +def lnprob_cz(pars, zprior, prior, zliml, zlimu, args, kwargs): + ''' + ''' + resid = residual_z(pars, *args, **kwargs) # i.e. (data - model) * wht + z = pars['z'] + s_z = 1 #pars['f_cz'] + resid *= 1/s_z + resid *= resid + + if False: + nzz = np.argmin(np.abs(zprior-z)) + + # For something unacceptable; + if nzz<0 or zprior[nzz]zlimu or prior[nzz]<=0: + return -np.inf + else: + respr = np.log(prior[nzz]) + resid += np.log(2 * np.pi * s_z**2) + return -0.5 * np.sum(resid) + respr + else: + resid += np.log(2 * np.pi * s_z**2) + return -0.5 * np.sum(resid) + + +def residual_z(pars, xm_tmp, fm_tmp, xobs, fobs, eobs, NR, NRbb_lim=10000, include_photometry=True, f_line_check=False): + ''' + ''' + vals = pars.valuesdict() + z = vals['z'] + Cz0s = vals['Cz0'] + Cz1s = vals['Cz1'] + Cz2s = vals['Cz2'] + + xm_s = xm_tmp * (1+z) + fint = interpolate.interp1d(xm_s, fm_tmp, kind='nearest', fill_value="extrapolate") + # @@@ Maybe for future consideration. + # fint = interpolate.interp1d(xm_s, fm_tmp, kind='linear', fill_value="extrapolate") + fm_s = fint(xobs) + + con0 = (NR<1000) + fy0 = fobs[con0] * Cz0s + ey0 = eobs[con0] * Cz0s + con1 = (NR>=1000) & (NR<2000) + fy1 = fobs[con1] * Cz1s + ey1 = eobs[con1] * Cz1s + con2 = (NR>=2000) & (NR=NRbb_lim) # BB + fy_bb = fobs[con_bb] + ey_bb = eobs[con_bb] + + fy01 = np.append(fy0,fy1) + ey01 = np.append(ey0,ey1) + fy02 = np.append(fy01,fy2) + ey02 = np.append(ey01,ey2) + + if include_photometry and len(fy_bb)>0: + # print('Including bb photometry in the redshift fit...') + fcon = np.append(fy02,fy_bb) + eycon = np.append(ey02,ey_bb) + wht = 1./np.square(eycon) + else: + # print('Not including bb photometry in the redshift fit...') + fcon = fy02 + eycon = ey02 + wht = 1./np.square(eycon) + + if f_line_check: + try: + wht2, ypoly = check_line_cz_man(fcon, xobs, wht, fm_s, z) + except: + wht2 = wht + else: + wht2 = wht + + if fobs is None: + print('Data is none') + return fm_s + else: + return (fm_s - fcon) * np.sqrt(wht2) # i.e. residual/sigma + def check_redshift(fobs, eobs, xobs, fm_tmp, xm_tmp, zbest, zprior, prior, NR, zliml, zlimu, \ - nmc_cz=100, nwalk_cz=10, nthin=5, f_line_check=False, f_vary=True, NRbb_lim=10000): + nmc_cz=100, nwalk_cz=10, nthin=5, f_line_check=False, f_vary=True, NRbb_lim=10000, include_photometry=True): ''' Purpose ------- @@ -39,11 +129,10 @@ def check_redshift(fobs, eobs, xobs, fm_tmp, xm_tmp, zbest, zprior, prior, NR, z fitc_cz : ''' - - from .function import check_line_cz_man - import numpy as np - from lmfit import Model, Parameters, minimize, fit_report, Minimizer - import scipy.interpolate as interpolate + if zliml == None or zlimu == None: + print('z range is not set for the z-fit function. Exiting.') + print('Specify `ZMCMIN` and `ZMCMIN` in your input file.') + sys.exit() fit_par_cz = Parameters() fit_par_cz.add('z', value=zbest, min=zliml, max=zlimu, vary=f_vary) @@ -51,77 +140,11 @@ def check_redshift(fobs, eobs, xobs, fm_tmp, xm_tmp, zbest, zprior, prior, NR, z fit_par_cz.add('Cz1', value=1, min=0.5, max=1.5) fit_par_cz.add('Cz2', value=1, min=0.5, max=1.5) - ############################## - def residual_z(pars): - vals = pars.valuesdict() - z = vals['z'] - Cz0s = vals['Cz0'] - Cz1s = vals['Cz1'] - Cz2s = vals['Cz2'] - - xm_s = xm_tmp * (1+z) - fint = interpolate.interp1d(xm_s, fm_tmp, kind='nearest', fill_value="extrapolate") - fm_s = fint(xobs) - - con0 = (NR<1000) - fy0 = fobs[con0] * Cz0s - ey0 = eobs[con0] * Cz0s - con1 = (NR>=1000) & (NR<2000) - fy1 = fobs[con1] * Cz1s - ey1 = eobs[con1] * Cz1s - con2 = (NR>=2000) & (NR=NRbb_lim) # BB - fy_bb = fobs[con_bb] - ey_bb = eobs[con_bb] - - fy01 = np.append(fy0,fy1) - ey01 = np.append(ey0,ey1) - fy02 = np.append(fy01,fy2) - ey02 = np.append(ey01,ey2) - - fcon = np.append(fy02,fy_bb) - eycon = np.append(ey02,ey_bb) - wht = 1./np.square(eycon) - - if f_line_check: - try: - wht2, ypoly = check_line_cz_man(fcon, xobs, wht, fm_s, z) - except: - wht2 = wht - else: - wht2 = wht - - if fobs is None: - print('Data is none') - return fm_s - else: - return (fm_s - fcon) * np.sqrt(wht2) # i.e. residual/sigma - - ############################### - def lnprob_cz(pars): - ''' - ''' - resid = residual_z(pars) # i.e. (data - model) * wht - z = pars['z'] - s_z = 1 #pars['f_cz'] - resid *= 1/s_z - resid *= resid - - nzz = np.argmin(np.abs(zprior-z)) - - # For something unacceptable; - if nzz<0 or zprior[nzz]zlimu or prior[nzz]<=0: - return -np.inf - else: - respr = np.log(prior[nzz]) - resid += np.log(2 * np.pi * s_z**2) - return -0.5 * np.sum(resid) + respr - ################################# - # Get Best fit - out_cz = minimize(residual_z, fit_par_cz, method='nelder') + args_res = (xm_tmp, fm_tmp, xobs, fobs, eobs, NR) + kwargs_res = {'include_photometry':include_photometry, 'f_line_check':f_line_check, 'NRbb_lim':NRbb_lim} + + out_cz = minimize(residual_z, fit_par_cz, args=args_res, method='nelder') keys = fit_report(out_cz).split('\n') for key in keys: if key[4:7] == 'chi': @@ -132,7 +155,38 @@ def lnprob_cz(pars): rcsq = float(skey[7]) fitc_cz = [csq, rcsq] # Chi2, Reduced-chi2 - mini_cz = Minimizer(lnprob_cz, out_cz.params) + + mini_cz = Minimizer(lnprob_cz, out_cz.params, (zprior, prior, zliml, zlimu, args_res, kwargs_res)) res_cz = mini_cz.emcee(burn=int(nmc_cz/2), steps=nmc_cz, thin=nthin, nwalkers=nwalk_cz, params=out_cz.params, is_weighted=True) return res_cz, fitc_cz + + +def get_chi2(zz_prob, fy_cz, ey_cz, x_cz, fm_tmp, xm_tmp, file_zprob, rms_lim=1e4): + ''' + zz_prob : float array + redshift array for fit. + fy_cz, ey_cz, x_cz : + observed values. + fm_tmp, xm_tmp : + template. + file_zprob : str + output file + ''' + mask = (ey_cz