From 45234e9da5a70b2efd9da78d04ca8741737d9ed2 Mon Sep 17 00:00:00 2001
From: Takahiro Morishita <tmorishita@stsci.edu>
Date: Tue, 25 Oct 2022 16:15:02 -0700
Subject: [PATCH 01/15] bug fixed

---
 gsf/fitting.py      | 30 ++++++++++++++++++++----------
 gsf/maketmp_filt.py | 24 +++++++++++++++++++-----
 gsf/zfit.py         | 25 ++++++++++++++++++++-----
 3 files changed, 59 insertions(+), 20 deletions(-)

diff --git a/gsf/fitting.py b/gsf/fitting.py
index bee457c..32f846b 100644
--- a/gsf/fitting.py
+++ b/gsf/fitting.py
@@ -183,12 +183,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
@@ -878,7 +881,7 @@ 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):
         '''
         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.
@@ -953,12 +956,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(0,zmax,delzz)
             if priors != None:
                 zprob = priors['z']
                 cprob = priors['chi2']
@@ -972,7 +975,7 @@ def fit_redshift(self, xm_tmp, fm_tmp, delzz=0.01, ezmin=0.01, zliml=0.01,
                     #prior_s /= np.sum(prior_s)
 
             else:
-                zz_prob = np.arange(0,13,delzz)
+                zz_prob = np.arange(0,zmax,delzz)
                 prior_s = zz_prob * 0 + 1.
                 prior_s /= np.sum(prior_s)
 
@@ -989,7 +992,10 @@ 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
@@ -1032,7 +1038,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.]
@@ -1067,7 +1073,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 +1099,17 @@ 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:
diff --git a/gsf/maketmp_filt.py b/gsf/maketmp_filt.py
index a24960f..e0d1b26 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')
diff --git a/gsf/zfit.py b/gsf/zfit.py
index 7086d32..a9a4539 100644
--- a/gsf/zfit.py
+++ b/gsf/zfit.py
@@ -1,6 +1,6 @@
 
 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
     -------
@@ -40,10 +40,16 @@ def check_redshift(fobs, eobs, xobs, fm_tmp, xm_tmp, zbest, zprior, prior, NR, z
 
     '''
 
-    from .function import check_line_cz_man
     import numpy as np
+    import sys
     from lmfit import Model, Parameters, minimize, fit_report, Minimizer
     import scipy.interpolate as interpolate
+    from .function import check_line_cz_man
+
+    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)
@@ -53,6 +59,8 @@ def check_redshift(fobs, eobs, xobs, fm_tmp, xm_tmp, zbest, zprior, prior, NR, z
 
     ##############################
     def residual_z(pars):
+        '''
+        '''
         vals = pars.valuesdict()
         z = vals['z']
         Cz0s = vals['Cz0']
@@ -81,9 +89,16 @@ def residual_z(pars):
         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 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:

From d7d99bf674164b5d672fe46b8c63795aa5889be1 Mon Sep 17 00:00:00 2001
From: Takahiro Morishita <takahiro@ipac.caltech.edu>
Date: Tue, 25 Oct 2022 17:17:49 -0700
Subject: [PATCH 02/15] scale

---
 gsf/fitting.py  |  7 +++----
 gsf/plot_sed.py | 30 +++++++++++++++++++++++-------
 2 files changed, 26 insertions(+), 11 deletions(-)

diff --git a/gsf/fitting.py b/gsf/fitting.py
index 32f846b..20d75ed 100644
--- a/gsf/fitting.py
+++ b/gsf/fitting.py
@@ -992,7 +992,7 @@ 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 = data_model[data_model[:, 0].argsort()]            
             #data_model_sort = np.sort(data_model, axis=0) # This does not work!!
 
@@ -1074,7 +1074,7 @@ def fit_redshift(self, xm_tmp, fm_tmp, delzz=0.01, ezmin=0.01, zliml=0.01,
             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 = data_model_new[data_model_new[:, 0].argsort()]            
+            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.)
@@ -1102,14 +1102,13 @@ def fit_redshift(self, xm_tmp, fm_tmp, delzz=0.01, ezmin=0.01, zliml=0.01,
             #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 = model_spec[model_spec[:, 0].argsort()]            
+                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:
diff --git a/gsf/plot_sed.py b/gsf/plot_sed.py
index 8d28f98..6d0498d 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
@@ -297,6 +297,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))) / 10
+    d = MB.d * scale
+
     #######################################
     # D.Kelson like Box for BB photometry
     #######################################
@@ -1198,12 +1204,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 +1408,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 +1619,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
     #######################################
@@ -2621,7 +2632,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 +2772,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 = (NR<NRbb_lim)
     sigma = 1.

From ee688f4540ced86fdd2a75df9bd80a2cef58e800 Mon Sep 17 00:00:00 2001
From: Takahiro Morishita <tmorishita@stsci.edu>
Date: Tue, 25 Oct 2022 17:33:22 -0700
Subject: [PATCH 03/15] Update plot_sed.py

---
 gsf/plot_sed.py | 16 ++++++++--------
 1 file changed, 8 insertions(+), 8 deletions(-)

diff --git a/gsf/plot_sed.py b/gsf/plot_sed.py
index 6d0498d..969f702 100644
--- a/gsf/plot_sed.py
+++ b/gsf/plot_sed.py
@@ -1061,10 +1061,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)
@@ -2322,10 +2322,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)

From 47135aece7b69d821da1a8e84f32a68f1291dc3d Mon Sep 17 00:00:00 2001
From: Takahiro Morishita <tmorishita@stsci.edu>
Date: Tue, 25 Oct 2022 22:22:41 -0700
Subject: [PATCH 04/15] zfit in progress

---
 gsf/fitting.py | 68 ++++++++++++++++++++++++++++++++++++++++----------
 gsf/zfit.py    | 43 +++++++++++++++++++++++++------
 2 files changed, 91 insertions(+), 20 deletions(-)

diff --git a/gsf/fitting.py b/gsf/fitting.py
index 20d75ed..9731768 100644
--- a/gsf/fitting.py
+++ b/gsf/fitting.py
@@ -21,7 +21,7 @@
 
 # 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
@@ -881,13 +881,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, zmax=20):
+        f_norm=True, f_lambda=False, zmax=20, fit_photometry=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
@@ -926,7 +932,11 @@ def fit_redshift(self, xm_tmp, fm_tmp, delzz=0.01, ezmin=0.01, zliml=0.01,
         sn = self.dict['fy'] / self.dict['ey']
 
         # Only spec data?
-        con_cz = (self.dict['NR']<NRbb_lim) & (sn>snlim)
+        if fit_photometry:
+            con_cz = (sn>snlim)
+        else:
+            con_cz = (self.dict['NR']<NRbb_lim) & (sn>snlim)
+
         if len(self.dict['fy'][con_cz])==0:
             if f_bb_zfit:
                 con_cz = (sn>snlim)
@@ -961,7 +971,8 @@ def fit_redshift(self, xm_tmp, fm_tmp, delzz=0.01, ezmin=0.01, zliml=0.01,
             prior_s = np.exp(-0.5 * cprob_s)
             prior_s /= np.sum(prior_s)
         else:
-            zz_prob = np.arange(0,zmax,delzz)
+            #zz_prob = np.arange(0,zmax,delzz)
+            zz_prob = np.arange(zliml,zlimu,delzz)
             if priors != None:
                 zprob = priors['z']
                 cprob = priors['chi2']
@@ -975,7 +986,7 @@ def fit_redshift(self, xm_tmp, fm_tmp, delzz=0.01, ezmin=0.01, zliml=0.01,
                     #prior_s /= np.sum(prior_s)
 
             else:
-                zz_prob = np.arange(0,zmax,delzz)
+                # zz_prob = np.arange(0,zmax,delzz)
                 prior_s = zz_prob * 0 + 1.
                 prior_s /= np.sum(prior_s)
 
@@ -1000,6 +1011,11 @@ def fit_redshift(self, xm_tmp, fm_tmp, delzz=0.01, ezmin=0.01, zliml=0.01,
             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:
+                self.file_zprob = self.DIR_OUT + 'zprob_' + self.ID + '.txt'
+                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('############################')
@@ -1168,6 +1184,8 @@ def fit_redshift(self, xm_tmp, fm_tmp, delzz=0.01, ezmin=0.01, zliml=0.01,
 
     def get_zdist(self, f_interact=False, f_ascii=True):
         '''
+        Purpose
+        -------
         Saves a plot of z-distribution.
 
         Parameters
@@ -1189,11 +1207,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]
@@ -1214,20 +1232,45 @@ 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
 
+        # 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):
         '''
@@ -1641,8 +1684,7 @@ 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
diff --git a/gsf/zfit.py b/gsf/zfit.py
index a9a4539..afed412 100644
--- a/gsf/zfit.py
+++ b/gsf/zfit.py
@@ -1,3 +1,9 @@
+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 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, include_photometry=True):
@@ -39,13 +45,6 @@ def check_redshift(fobs, eobs, xobs, fm_tmp, xm_tmp, zbest, zprior, prior, NR, z
     fitc_cz :
 
     '''
-
-    import numpy as np
-    import sys
-    from lmfit import Model, Parameters, minimize, fit_report, Minimizer
-    import scipy.interpolate as interpolate
-    from .function import check_line_cz_man
-
     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.')
@@ -151,3 +150,33 @@ def lnprob_cz(pars):
     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<rms_lim)
+    prob_cz = np.zeros(len(zz_prob), float)
+    fw = open(file_zprob, 'w')
+    fw.write('# z p(z)\n')
+    for ii in range(len(zz_prob)):
+        z = zz_prob[ii]
+        xm_s = xm_tmp * (1+z)
+        fint = interpolate.interp1d(xm_s, fm_tmp, kind='nearest', fill_value="extrapolate")
+        fm_s = fint(x_cz)
+
+        wht = 1./np.square(ey_cz)
+        # print(np.nansum((fm_s - fy_cz)**2 * (wht)))
+        lnprob_cz = -0.5 * np.nansum( np.square((fm_s - fy_cz)[mask] * np.sqrt(wht[mask])) ) # i.e. (residual/sigma)^2
+        prob_cz[ii] = np.exp(lnprob_cz)
+        fw.write('%.3f %.3e\n'%(z,prob_cz[ii]))
+        
+    fw.close()

From babf571c7ac612f052dd427bc3736c3319eb0ba6 Mon Sep 17 00:00:00 2001
From: Takahiro Morishita <takahiro@ipac.caltech.edu>
Date: Wed, 26 Oct 2022 11:36:21 -0700
Subject: [PATCH 05/15] cleaned the code

---
 gsf/fitting.py |  70 +++++++++++----------
 gsf/zfit.py    | 164 +++++++++++++++++++++++++------------------------
 2 files changed, 120 insertions(+), 114 deletions(-)

diff --git a/gsf/fitting.py b/gsf/fitting.py
index 9731768..41f6cd9 100644
--- a/gsf/fitting.py
+++ b/gsf/fitting.py
@@ -16,6 +16,7 @@
 import timeit
 from scipy import stats
 from scipy.stats import norm
+import scipy.interpolate as interpolate
 from astropy.io import fits,ascii
 import corner
 
@@ -802,7 +803,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')
@@ -881,7 +881,7 @@ 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, zmax=20, fit_photometry=False):
+        f_norm=True, f_lambda=False, zmax=20, fit_photometry=False, f_exclude_negative=False):
         '''
         Purpose
         -------
@@ -910,21 +910,20 @@ 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
 
@@ -933,16 +932,19 @@ def fit_redshift(self, xm_tmp, fm_tmp, delzz=0.01, ezmin=0.01, zliml=0.01,
 
         # Only spec data?
         if fit_photometry:
-            con_cz = (sn>snlim)
+            con_cz = ()#(sn>snlim)
         else:
-            con_cz = (self.dict['NR']<NRbb_lim) & (sn>snlim)
+            con_cz = (self.dict['NR']<NRbb_lim) #& (sn>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
@@ -1013,7 +1015,6 @@ def fit_redshift(self, xm_tmp, fm_tmp, delzz=0.01, ezmin=0.01, zliml=0.01,
 
             # Write prob distribution;
             if True:
-                self.file_zprob = self.DIR_OUT + 'zprob_' + self.ID + '.txt'
                 get_chi2(zz_prob, fy_cz, ey_cz, x_cz, fm_tmp, xm_tmp/(1+self.zgal), self.file_zprob)
 
             print('############################')
@@ -1246,30 +1247,30 @@ def get_zdist(self, f_interact=False, f_ascii=True):
             print('z-distribution figure is not generated.')
             pass
 
-        # 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()
 
+        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):
@@ -1598,7 +1599,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):
         '''
         Main module of this script.
 
@@ -1767,7 +1769,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)
diff --git a/gsf/zfit.py b/gsf/zfit.py
index afed412..6cd0b99 100644
--- a/gsf/zfit.py
+++ b/gsf/zfit.py
@@ -5,6 +5,84 @@
 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
+    
+    nzz = np.argmin(np.abs(zprior-z))
+
+    # For something unacceptable;
+    if nzz<0 or zprior[nzz]<zliml 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
+
+
+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")
+    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)
+    fy2  = fobs[con2] * Cz2s
+    ey2  = eobs[con2] * Cz2s
+    con_bb = (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, include_photometry=True):
     '''
@@ -56,86 +134,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)
-        fy2  = fobs[con2] * Cz2s
-        ey2  = eobs[con2] * Cz2s
-        con_bb = (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 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]<zliml 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':
@@ -146,7 +149,8 @@ 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

From 4414958e29bd9cf112146127c79e0816cc63b0bb Mon Sep 17 00:00:00 2001
From: Takahiro Morishita <takahiro@ipac.caltech.edu>
Date: Wed, 26 Oct 2022 15:38:17 -0700
Subject: [PATCH 06/15] logufix

---
 gsf/fitting.py    | 2 +-
 gsf/maketmp_z0.py | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/gsf/fitting.py b/gsf/fitting.py
index 41f6cd9..284be9c 100644
--- a/gsf/fitting.py
+++ b/gsf/fitting.py
@@ -231,7 +231,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:
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.

From f93728f2456dd22f443c24a0801c587956672e9f Mon Sep 17 00:00:00 2001
From: Takahiro Morishita <takahiro@ipac.caltech.edu>
Date: Wed, 26 Oct 2022 18:17:21 -0700
Subject: [PATCH 07/15] zfit improving

---
 gsf/fitting.py | 13 ++++++++++---
 gsf/writing.py |  2 +-
 gsf/zfit.py    | 18 +++++++++++-------
 3 files changed, 22 insertions(+), 11 deletions(-)

diff --git a/gsf/fitting.py b/gsf/fitting.py
index 284be9c..8af78c1 100644
--- a/gsf/fitting.py
+++ b/gsf/fitting.py
@@ -1021,7 +1021,7 @@ def fit_redshift(self, xm_tmp, fm_tmp, delzz=0.01, ezmin=0.01, zliml=0.01,
             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])
@@ -1169,6 +1169,11 @@ def fit_redshift(self, xm_tmp, fm_tmp, delzz=0.01, ezmin=0.01, zliml=0.01,
         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
@@ -1693,8 +1698,6 @@ def main(self, cornerplot:bool=True, specplot=1, sigz=1.0, ezmin=0.01, ferr=0,
         #################################################
         if flag_z == 'y' or flag_z == '':
 
-            self.get_zdist()
-
             #######################
             # Add parameters;
             #######################
@@ -2057,8 +2060,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))
+
+            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/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 6cd0b99..953f2e2 100644
--- a/gsf/zfit.py
+++ b/gsf/zfit.py
@@ -14,16 +14,20 @@ def lnprob_cz(pars, zprior, prior, zliml, zlimu, args, kwargs):
     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]<zliml or zprior[nzz]>zlimu or prior[nzz]<=0:
-        return -np.inf
+    if False:
+        nzz = np.argmin(np.abs(zprior-z))
+
+        # For something unacceptable;
+        if nzz<0 or zprior[nzz]<zliml 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:
-        respr = np.log(prior[nzz])
         resid += np.log(2 * np.pi * s_z**2)
-        return -0.5 * np.sum(resid) + respr
+        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):

From 5dd26dba5081dca8498043685e5a608263ee71cb Mon Sep 17 00:00:00 2001
From: Takahiro Morishita <takahiro@ipac.caltech.edu>
Date: Wed, 26 Oct 2022 21:03:48 -0700
Subject: [PATCH 08/15] redshift fit

---
 gsf/fitting.py      | 23 ++++++++++++++++-------
 gsf/maketmp_filt.py |  2 +-
 gsf/plot_sed.py     | 26 ++++++++++++--------------
 gsf/zfit.py         |  2 ++
 4 files changed, 31 insertions(+), 22 deletions(-)

diff --git a/gsf/fitting.py b/gsf/fitting.py
index 8af78c1..83ca9aa 100644
--- a/gsf/fitting.py
+++ b/gsf/fitting.py
@@ -94,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:
@@ -409,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:
@@ -421,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:
@@ -441,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
@@ -1042,6 +1047,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.')
@@ -1164,8 +1170,8 @@ 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'
 
@@ -1184,6 +1190,7 @@ 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
 
@@ -1447,10 +1454,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)
@@ -1668,6 +1675,8 @@ def main(self, cornerplot:bool=True, specplot=1, sigz=1.0, ezmin=0.01, ferr=0,
             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('#####################################')
@@ -2061,8 +2070,8 @@ 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()
 
diff --git a/gsf/maketmp_filt.py b/gsf/maketmp_filt.py
index e0d1b26..d23ab07 100644
--- a/gsf/maketmp_filt.py
+++ b/gsf/maketmp_filt.py
@@ -631,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/plot_sed.py b/gsf/plot_sed.py
index 969f702..39e16c8 100644
--- a/gsf/plot_sed.py
+++ b/gsf/plot_sed.py
@@ -784,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:
@@ -1018,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
diff --git a/gsf/zfit.py b/gsf/zfit.py
index 953f2e2..64271b8 100644
--- a/gsf/zfit.py
+++ b/gsf/zfit.py
@@ -41,6 +41,8 @@ def residual_z(pars, xm_tmp, fm_tmp, xobs, fobs, eobs, NR, NRbb_lim=10000, inclu
 
     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)

From c16fe5f1acaf3307db69e3541e64cf29ebef5237 Mon Sep 17 00:00:00 2001
From: Takahiro Morishita <takahiro@ipac.caltech.edu>
Date: Thu, 27 Oct 2022 07:26:33 -0700
Subject: [PATCH 09/15] spec only option for z search

---
 gsf/fitting.py            |  6 ++++--
 gsf/function.py           |  8 ++++++--
 gsf/posterior_flexible.py | 10 +++++++++-
 3 files changed, 19 insertions(+), 5 deletions(-)

diff --git a/gsf/fitting.py b/gsf/fitting.py
index 83ca9aa..f9122ed 100644
--- a/gsf/fitting.py
+++ b/gsf/fitting.py
@@ -1666,9 +1666,11 @@ 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
diff --git a/gsf/function.py b/gsf/function.py
index d02a85f..a7bb33f 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
diff --git a/gsf/posterior_flexible.py b/gsf/posterior_flexible.py
index e2f5948..2116b99 100644
--- a/gsf/posterior_flexible.py
+++ b/gsf/posterior_flexible.py
@@ -19,7 +19,8 @@ def __init__(self, mainbody):
         self.Na = len(self.mb.age)
 
 
-    def residual(self, pars, fy:float, ey:float, wht:float, f_fir:bool, out:bool=False, f_val:bool=False, f_penlize:bool=True):
+    def residual(self, pars, fy:float, ey:float, wht:float, f_fir:bool, out:bool=False, 
+        f_val:bool=False, f_penlize:bool=True, f_only_spec:bool=False):
         '''
         Parameters
         ----------
@@ -62,6 +63,13 @@ def residual(self, pars, fy:float, ey:float, wht:float, f_fir:bool, out:bool=Fal
         else:
             logf = -np.inf # temporary... (if f is param, then take from vals dictionary.)
 
+        if f_only_spec:
+            contmp = (self.mb.dict['NR']<10000)
+            fy = fy[contmp]
+            ey = ey[contmp]
+            wht = wht[contmp]
+            model = model[contmp]
+
         sig = wht[:] * 0
         con_res = (wht>0)
         con_res_r = (wht==0)

From dbf49142f9caf2135e1fbd5997173692a91a0141 Mon Sep 17 00:00:00 2001
From: Takahiro Morishita <takahiro@ipac.caltech.edu>
Date: Sat, 29 Oct 2022 15:16:43 -0700
Subject: [PATCH 10/15] bug fixed

---
 gsf/fitting.py | 20 +++++++++++---------
 1 file changed, 11 insertions(+), 9 deletions(-)

diff --git a/gsf/fitting.py b/gsf/fitting.py
index f9122ed..3cf6f36 100644
--- a/gsf/fitting.py
+++ b/gsf/fitting.py
@@ -1047,7 +1047,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]
+            #self.fitc_cz = fitc_cz[1]
 
         else:
             print('fzvis is set to False. z fit not happening.')
@@ -1069,6 +1069,8 @@ 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
+            #self.fitc_cz = z_cz[1]
+            fitc_cz = [99,99]
 
             # If this label is being used, it means that the fit is failed.
             fit_label = 'Current model'
@@ -1293,10 +1295,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
@@ -1612,7 +1614,7 @@ 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_chind:bool=True):
+            f_chind:bool=True, ncpu:int=0):
         '''
         Main module of this script.
 
@@ -1727,8 +1729,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:
@@ -1752,7 +1755,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 ##')

From 2958cdf4ca0b3f2c275f289f40353395154dd298 Mon Sep 17 00:00:00 2001
From: Takahiro Morishita <takahiro@ipac.caltech.edu>
Date: Sat, 29 Oct 2022 22:25:13 -0700
Subject: [PATCH 11/15] Update plot_sfh.py

---
 gsf/plot_sfh.py | 32 ++++++++++++++++++++++++--------
 1 file changed, 24 insertions(+), 8 deletions(-)

diff --git a/gsf/plot_sfh.py b/gsf/plot_sfh.py
index cc3369f..2fcd1ff 100644
--- a/gsf/plot_sfh.py
+++ b/gsf/plot_sfh.py
@@ -538,9 +538,12 @@ def plot_sfh(MB, f_comp=0, flim=0.01, lsfrl=-3, mmax=1000, Txmin=0.08, Txmax=4,
     elif zbes<6:
         zred  = [zbes, 5, 6, 9]
         zredl = ['$z_\mathrm{obs.}$', 5, 6, 9]
-    else:
+    elif zbes<12:
         zred  = [zbes, 12]
         zredl = ['$z_\mathrm{obs.}$', 12]
+    else:
+        zred  = [zbes, 20]
+        zredl = ['$z_\mathrm{obs.}$', 20]
 
     Tzz = np.zeros(len(zred), dtype='float')
     for zz in range(len(zred)):
@@ -652,7 +655,7 @@ def plot_sfh(MB, f_comp=0, flim=0.01, lsfrl=-3, mmax=1000, Txmin=0.08, Txmax=4,
     col50 = fits.Column(name='Z84', format='E', unit='logZsun', array=ZCp[:,2])
     col02.append(col50)
     
-    colms  = fits.ColDefs(col02)
+    colms = fits.ColDefs(col02)
     dathdu = fits.BinTableHDU.from_columns(colms)
     hdu = fits.HDUList([prihdu, dathdu])
     file_sfh = MB.DIR_OUT + 'SFH_' + ID + '.fits'
@@ -685,17 +688,27 @@ def plot_sfh(MB, f_comp=0, flim=0.01, lsfrl=-3, mmax=1000, Txmin=0.08, Txmax=4,
     ax2.set_yticklabels(np.arange(y2min, y2max, dely2), minor=False)
 
     ax2.yaxis.set_major_formatter(FormatStrFormatter('%.1f'))
-    y3min, y3max = np.min(Z), np.max(Z)
 
     if not skip_zhist:
+
+        delZl = 0.4 # only for label. not delZ of Zall.
+        delZZl = delZl / 4
+
+        y3min, y3max = np.min(Z), np.max(Z)+delZZl
         ax4.set_xlim(Txmin, Txmax)
-        ax4.set_ylim(y3min-0.05, y3max)
+        # There is a bug that lower size does not match with label.
+        # ax4.set_ylim(y3min-delZZl, y3max+delZZl)
+        ax4.set_ylim(y3min, y3max)
         ax4.set_xscale('log')
 
         if f_axis_force:
-            ax4.set_yticks([-0.8, -0.4, 0., 0.4])
-            ax4.set_yticklabels(['-0.8', '-0.4', '0', '0.4'])
-        
+            Zticks = np.arange(MB.Zall.min(), MB.Zall.max()+delZZl, delZl)
+            ax4.set_yticks(Zticks)
+            Zlabels = []
+            for Ztmp in Zticks:
+                Zlabels.append('%.1f'%(Ztmp))
+            #ax4.set_yticklabels(Zlabels)
+
         ax4.yaxis.labelpad = -2
         ax4.set_xlabel('$t_\mathrm{lookback}$/Gyr', fontsize=12)
         ax4.set_ylabel('$\log Z_*/Z_\odot$', fontsize=12)
@@ -1205,9 +1218,12 @@ def plot_sfh_tau(MB, f_comp=0, flim=0.01, lsfrl=-1, mmax=1000, Txmin=0.08, Txmax
     elif zbes<6:
         zred  = [zbes, 5, 6, 9]
         zredl = ['$z_\mathrm{obs.}$', 5, 6, 9]
-    else:
+    elif zbes<12:
         zred  = [zbes, 12]
         zredl = ['$z_\mathrm{obs.}$', 12]
+    else:
+        zred  = [zbes, 20]
+        zredl = ['$z_\mathrm{obs.}$', 20]
 
     Tzz = np.zeros(len(zred), dtype='float')
     for zz in range(len(zred)):

From 029de876eb1e5c16914f7d2646f14f40780c4120 Mon Sep 17 00:00:00 2001
From: Takahiro Morishita <takahiro@ipac.caltech.edu>
Date: Sun, 30 Oct 2022 12:37:05 -0700
Subject: [PATCH 12/15] Update plot_sed.py

---
 gsf/plot_sed.py | 12 ++++++------
 1 file changed, 6 insertions(+), 6 deletions(-)

diff --git a/gsf/plot_sed.py b/gsf/plot_sed.py
index 39e16c8..d0e867a 100644
--- a/gsf/plot_sed.py
+++ b/gsf/plot_sed.py
@@ -289,12 +289,12 @@ 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;
@@ -501,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

From d311745829972633790f76c9fb19c9d9de6146b6 Mon Sep 17 00:00:00 2001
From: Takahiro Morishita <takahiro@ipac.caltech.edu>
Date: Sat, 12 Nov 2022 15:39:29 -0800
Subject: [PATCH 13/15] Update function.py

---
 gsf/function.py | 12 +++++++++---
 1 file changed, 9 insertions(+), 3 deletions(-)

diff --git a/gsf/function.py b/gsf/function.py
index a7bb33f..83212cd 100644
--- a/gsf/function.py
+++ b/gsf/function.py
@@ -554,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
 

From ab9d84c9f18f52029752ef23bdbdb771251b0b86 Mon Sep 17 00:00:00 2001
From: Takahiro Morishita <takahiro@ipac.caltech.edu>
Date: Sat, 19 Nov 2022 15:10:28 -0500
Subject: [PATCH 14/15] docs

---
 docs/parameters.rst   | 5 ++++-
 gsf/function_class.py | 4 ++--
 2 files changed, 6 insertions(+), 3 deletions(-)

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/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)

From 9cf5fef0f37149512f8855ec6945074bcc421d94 Mon Sep 17 00:00:00 2001
From: Takahiro Morishita <takahiro@ipac.caltech.edu>
Date: Tue, 22 Nov 2022 16:30:08 -0800
Subject: [PATCH 15/15] Update fitting.py

---
 gsf/fitting.py | 15 ++++++++-------
 1 file changed, 8 insertions(+), 7 deletions(-)

diff --git a/gsf/fitting.py b/gsf/fitting.py
index 3cf6f36..1d34656 100644
--- a/gsf/fitting.py
+++ b/gsf/fitting.py
@@ -929,8 +929,14 @@ def fit_redshift(self, xm_tmp, fm_tmp, delzz=0.01, ezmin=0.01, zliml=0.01,
         self.nmc_cz = int(self.inputs['NMCZ'])
 
         # For z prior.
-        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']
@@ -978,7 +984,6 @@ def fit_redshift(self, xm_tmp, fm_tmp, delzz=0.01, ezmin=0.01, zliml=0.01,
             prior_s = np.exp(-0.5 * cprob_s)
             prior_s /= np.sum(prior_s)
         else:
-            #zz_prob = np.arange(0,zmax,delzz)
             zz_prob = np.arange(zliml,zlimu,delzz)
             if priors != None:
                 zprob = priors['z']
@@ -990,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,zmax,delzz)
                 prior_s = zz_prob * 0 + 1.
                 prior_s /= np.sum(prior_s)
 
@@ -1069,7 +1071,6 @@ 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
-            #self.fitc_cz = z_cz[1]
             fitc_cz = [99,99]
 
             # If this label is being used, it means that the fit is failed.