diff --git a/gsf/fitting.py b/gsf/fitting.py
index d19bd7e..88773e7 100644
--- a/gsf/fitting.py
+++ b/gsf/fitting.py
@@ -992,15 +992,17 @@ def fit_redshift(self, xm_tmp, fm_tmp, delzz=0.01, ezmin=0.01, #zliml=0.01, zlim
             # Observed data;
             # spec;
             if self.has_spectrum:
-                con = (self.dict['ey']<1000) & (self.dict['NR']<self.NRbb_lim)
+                ey_max = 1000
+                con = (self.dict['ey']<ey_max) & (self.dict['NR']<self.NRbb_lim) & (self.dict['ey']>=0)
                 ax1.errorbar(self.dict['x'][con], self.dict['fy'][con], yerr=self.dict['ey'][con], color='gray', capsize=0, linewidth=0.5, linestyle='', zorder=4)
                 ax1.plot(self.dict['x'][con], self.dict['fy'][con], '.r', linestyle='', linewidth=0.5, label='Observed spectrum', zorder=4)
             # bb;
-            con = (self.dict['NR']>=self.NRbb_lim)
-            ax1.errorbar(self.dict['x'][con], self.dict['fy'][con], yerr=self.dict['ey'][con], ms=15, marker='None', 
-                color='orange', capsize=0, linewidth=0.5, ls='None', label='', zorder=4)
-            ax1.scatter(self.dict['x'][con], self.dict['fy'][con], s=100, marker='o', 
-                color='orange', edgecolor='k', label='Observed photometry', zorder=4)
+            if include_photometry:
+                con = (self.dict['NR']>=self.NRbb_lim) & (self.dict['ey']>=0)
+                ax1.errorbar(self.dict['x'][con], self.dict['fy'][con], yerr=self.dict['ey'][con], ms=15, marker='None', 
+                    color='orange', capsize=0, linewidth=0.5, ls='None', label='', zorder=4)
+                ax1.scatter(self.dict['x'][con], self.dict['fy'][con], s=100, marker='o', 
+                    color='orange', edgecolor='k', label='Observed photometry', zorder=4)
 
             # Write prob distribution;
             get_chi2(zz_prob, fy_cz, ey_cz, x_cz, fm_tmp, xm_tmp/(1+self.zgal), self.file_zprob)
@@ -1322,9 +1324,7 @@ def set_param(self):
         '''
         Set parameters.
         '''
-        print('##################')
-        print('Setting parameters')
-        print('##################\n')
+        self.logger.info('Setting parameters')
         agemax = self.cosmo.age(self.zgal).value 
         fit_params = Parameters()
         f_Alog = True
@@ -1482,9 +1482,7 @@ def check_mainbody(self):
     def prepare_class(self, add_fir=None):
         '''
         '''
-        print('#################')
-        print('Preparing library')
-        print('#################\n')
+        self.logger.info('Preparing library')
         # Load Spectral library;
         self.lib = self.fnc.open_spec_fits(fall=0)
         self.lib_all = self.fnc.open_spec_fits(fall=1, orig=True)
@@ -1702,21 +1700,11 @@ def main(self, cornerplot:bool=True, specplot=1, sigz=1.0, ezmin=0.01, ferr=0,
                     except:
                         out.params['Z0'].value = out_keep.params['Z0'].value
 
-
             ##############################
-            print('\n\n')
-            print('###############################')
-            print('Input redshift is adopted.')
-            print('Starting long journey in MCMC.')
-            print('###############################')
-            print('\n\n')
-
-            ################################
-            self.logger.info('Minimizer Defined')
-
-            print('########################')
-            print('### Starting sampling ##')
-            print('########################\n')
+            self.logger.info('Input redshift is adopted')
+            self.logger.info('Starting MCMC')
+            # self.logger.info('Minimizer Defined')
+            # self.logger.info('Starting sampling')
             start_mc = timeit.default_timer()
 
             # MCMC;
diff --git a/gsf/gsf.py b/gsf/gsf.py
index 7165460..5c54fac 100644
--- a/gsf/gsf.py
+++ b/gsf/gsf.py
@@ -87,7 +87,7 @@ def run_gsf_all(parfile, fplt, cornerplot=True, f_plot_chain=True, f_Alog=True,
     f_prior_sfh=False, norder_sfh_prior=3, 
     f_shuffle=False, amp_shuffle=1e-2, Zini=None, tau_lim=0.001,
     skip_sfh=False, f_fancyplot=False, skip_zhist=False, f_sfh_yaxis_force=True, tset_SFR_SED=0.1, 
-    nthin=1, delwave=1, f_plot_resid=False, scale=1e-19, f_plot_filter=True,
+    nthin=1, delwave=1, f_plot_resid=False, scale=None, f_plot_filter=True,
     mmax_param:int=1000
     ):
     '''
diff --git a/gsf/maketmp_filt.py b/gsf/maketmp_filt.py
index f11e012..317440b 100644
--- a/gsf/maketmp_filt.py
+++ b/gsf/maketmp_filt.py
@@ -141,12 +141,11 @@ def check_library(MB, af, nround=3):
     
     flag = True
     MB.logger.info('Checking the template library...')
-    print('Speficied - Template')
 
     # No. of age;
     if MB.SFH_FORM==-99:
         if len(af['ML']['ms_0']) != len(MB.age):
-            print('No of age pixels:', len(MB.age), len(af['ML']['ms_0']))
+            MB.logger.error('No of age pixels:', len(MB.age), len(af['ML']['ms_0']))
             flag = False
     else:
         flag = True
@@ -154,36 +153,39 @@ def check_library(MB, af, nround=3):
     # Matallicity:
     for aa in range(len(Zall)):
         if Zall[aa] != af['Z%d'%(aa)]:
-            print('Z:', Zall[aa], af['Z%d'%(aa)])
+            MB.logger.error('Z:', Zall[aa], af['Z%d'%(aa)])
             flag = False
 
     if MB.SFH_FORM==-99:
         # Age:
         for aa in range(len(MB.age)):
             if round(MB.age[aa],nround) != round(af['age%d'%(aa)],nround):
-                print('age:', MB.age[aa], af['age%d'%(aa)])
+                MB.logger.error('age:', MB.age[aa], af['age%d'%(aa)])
                 flag = False
         # Tau (e.g. ssp/csp):
         for aa in range(len(MB.tau0)):
             if round(MB.tau0[aa]) != round(af['tau0%d'%(aa)]):
-                print('tau0:', MB.tau0[aa], af['tau0%d'%(aa)])
+                MB.logger.error('tau0:', MB.tau0[aa], af['tau0%d'%(aa)])
                 flag = False
     else:
         # Age:
         for aa in range(len(MB.ageparam)):
             if round(MB.ageparam[aa]) != round(af['age%d'%(aa)]):
-                print('age:', MB.ageparam[aa], af['age%d'%(aa)])
+                MB.logger.error('age:', MB.ageparam[aa], af['age%d'%(aa)])
                 flag = False
         for aa in range(len(MB.tau)):
             if round(MB.tau[aa]) != round(af['tau%d'%(aa)]):
-                print('tau:', MB.tau[aa], af['tau%d'%(aa)])
+                MB.logger.error('tau:', MB.tau[aa], af['tau%d'%(aa)])
                 flag = False
 
     # IMF:
     if MB.nimf != af['nimf']:
-        print('nimf:', MB.nimf, af['nimf'])
+        MB.logger.error('nimf:', MB.nimf, af['nimf'])
         flag = False
 
+    if not flag:
+        MB.logger.error('# Specified - Template')
+
     return flag
 
 
@@ -512,64 +514,91 @@ def maketemp(MB, ebblim=1e10, lamliml=0., lamlimu=50000., ncolbb=10000,
 
     try:
         spec_files = [x.strip() for x in inputs['SPEC_FILE'].split(',')]
-        ninp0 = np.zeros(len(spec_files), dtype='int')
+    except:
+        spec_files = []
+        MB.logger.info('No spec file is provided.')
+        pass
+    ninp0 = np.zeros(len(spec_files), dtype='int')
 
-        # THIS PART IS JUST TO GET THE TOTAL ARRAY NUMBER;
-        for ff, spec_file in enumerate(spec_files):
-            try:
-                if spec_file.split('.')[-1] == 'asdf':
-                    id_asdf = int(spec_file.split('_')[2])
-                    fd0 = asdf.open(os.path.join(DIR_EXTR, spec_file))
-                    lm0tmp = fd0[id_asdf]['wavelength'].to(u.angstrom)
-                    ninp0[ff] = len(lm0tmp)
-                else:
-                    fd0 = np.loadtxt(os.path.join(DIR_EXTR, spec_file), comments='#')
-                    lm0tmp = fd0[:,0]
-                    ninp0[ff] = len(lm0tmp)
-            except Exception:
-                MB.logger.warning('File, %s, cannot be open.'%(os.path.join(DIR_EXTR, spec_file)))
-                pass
+    # THIS PART IS JUST TO GET THE TOTAL ARRAY NUMBER;
+    for ff, spec_file in enumerate(spec_files):
+        try:
+            if spec_file.split('.')[-1] == 'asdf':
+                id_asdf = int(spec_file.split('_')[2])
+                fd0 = asdf.open(os.path.join(DIR_EXTR, spec_file))
+                lm0tmp = fd0[id_asdf]['wavelength'].to(u.angstrom)
+                ninp0[ff] = len(lm0tmp)
+            elif spec_file.split('.')[-1] == 'fits':
+                fd0 = fits.open(os.path.join(DIR_EXTR, spec_file))[1].data
+                eobs0 = fd0['full_err']
+                spec_mask = (eobs0>0)
+                lm0tmp = fd0['wave'][spec_mask]
+                ninp0[ff] = len(lm0tmp)
+            else:
+                fd0 = np.loadtxt(os.path.join(DIR_EXTR, spec_file), comments='#')
+                lm0tmp = fd0[:,0]
+                ninp0[ff] = len(lm0tmp)
+        except Exception:
+            MB.logger.error('File, %s, cannot be open.'%(os.path.join(DIR_EXTR, spec_file)))
+            pass
 
-        # Then, Constructing arrays.
-        lm = np.zeros(np.sum(ninp0[:]),dtype=float)
-        fobs = np.zeros(np.sum(ninp0[:]),dtype=float)
-        eobs = np.zeros(np.sum(ninp0[:]),dtype=float)
-        fgrs = np.zeros(np.sum(ninp0[:]),dtype=int) # FLAG for each grism.
-        for ff, spec_file in enumerate(spec_files):
-            if True:#try:
-                if spec_file.split('.')[-1] == 'asdf':
-                    id_asdf = int(spec_file.split('_')[2])
-                    fd0 = asdf.open(os.path.join(DIR_EXTR, spec_file))
-                    lm0tmp = fd0[id_asdf]['wavelength'].to(u.angstrom).value
-                    fobs0 = fd0[id_asdf]['flux'].value
-                    eobs0 = np.sqrt(fd0[id_asdf]['fluxvar']).value
+    # Then, Constructing arrays.
+    lm = np.zeros(np.sum(ninp0[:]),dtype=float)
+    fobs = np.zeros(np.sum(ninp0[:]),dtype=float)
+    eobs = np.zeros(np.sum(ninp0[:]),dtype=float)
+    fgrs = np.zeros(np.sum(ninp0[:]),dtype=int) # FLAG for each grism.
+    for ff, spec_file in enumerate(spec_files):
+        try:
+            if spec_file.split('.')[-1] == 'asdf':
+                id_asdf = int(spec_file.split('_')[2])
+                fd0 = asdf.open(os.path.join(DIR_EXTR, spec_file))
+                lm0tmp = fd0[id_asdf]['wavelength'].to(u.angstrom).value
+                fobs0 = fd0[id_asdf]['flux'].value
+                eobs0 = np.sqrt(fd0[id_asdf]['fluxvar']).value
+            elif spec_file.split('.')[-1] == 'fits':
+                fd0 = fits.open(os.path.join(DIR_EXTR, spec_file))[1].data
+                eobs0 = fd0['full_err']
+                if True:
+                    spec_mask = (eobs0>0)
                 else:
-                    fd0 = np.loadtxt(os.path.join(DIR_EXTR, spec_file), comments='#')
-                    lm0tmp = fd0[:,0]
-                    fobs0 = fd0[:,1]
-                    eobs0 = fd0[:,2]
-
-                for ii1 in range(ninp0[ff]):
-                    if ff==0:
-                        ii = ii1
-                    else:
-                        ii = ii1 + np.sum(ninp0[:ff])
-                    fgrs[ii] = ff
-                    lm[ii] = lm0tmp[ii1]
-                    fobs[ii] = fobs0[ii1]
-                    eobs[ii] = eobs0[ii1]
-
-                MB.f_spec = True
-                data_meta['data_len'][ff] = len(lm0tmp)
-                data_meta['data_origin'] = np.append(data_meta['data_origin'], '%s'%spec_file)
-                data_meta['data_index'] = np.append(data_meta['data_index'], '%d'%ff)
+                    spec_mask = ()
+                lm0tmp = fd0['wave'][spec_mask]
+                if lm0tmp.max() < 10:
+                    MB.logger.warning('Wave column in the input spec file seems to be um. Scaling to AA.')
+                    lm0tmp *= 1e4
+                try:
+                    magzp_spec = float(inputs['MAGZP_SPEC'])
+                    magzp = float(inputs['MAGZP'])
+                    C_spec = 10**((magzp-magzp_spec)/(2.5))
+                except:
+                    MB.logger.info('`MAGZP_SPEC` not found. Assuming same as `MAGZP`')
+                    C_spec = 1
+                fobs0 = fd0['flux'][spec_mask] * C_spec
+                eobs0 = fd0['full_err'][spec_mask] * C_spec
+            else:
+                fd0 = np.loadtxt(os.path.join(DIR_EXTR, spec_file), comments='#')
+                lm0tmp = fd0[:,0]
+                fobs0 = fd0[:,1]
+                eobs0 = fd0[:,2]
 
-            else:#except Exception:
-                print('No spec data is registered.')
-                pass
-    except:
-        MB.logger.info('No spec file is provided.')
-        pass
+            for ii1 in range(ninp0[ff]):
+                if ff==0:
+                    ii = ii1
+                else:
+                    ii = ii1 + np.sum(ninp0[:ff])
+                fgrs[ii] = ff
+                lm[ii] = lm0tmp[ii1]
+                fobs[ii] = fobs0[ii1]
+                eobs[ii] = eobs0[ii1]
+
+            data_meta['data_len'][ff] = len(lm0tmp)
+            data_meta['data_origin'] = np.append(data_meta['data_origin'], '%s'%spec_file)
+            data_meta['data_index'] = np.append(data_meta['data_index'], '%d'%ff)
+            MB.f_spec = True
+
+        except Exception:
+            print('No spec data is registered.')
+            pass
 
     if ncolbb < np.sum(data_meta['data_len']):
         MB.logger.info('ncolbb is updated')
diff --git a/gsf/plot_sed.py b/gsf/plot_sed.py
index bf624fb..cdb55c5 100644
--- a/gsf/plot_sed.py
+++ b/gsf/plot_sed.py
@@ -33,10 +33,10 @@
 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, 
+def plot_sed(MB, flim=0.01, fil_path='./', scale=None, 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,
-    x1min=4000, return_figure=False):
+    x1min=4000, return_figure=False, lcb='#4682b4'):
     '''
     Parameters
     ----------
@@ -56,25 +56,7 @@ def plot_sed(MB, flim=0.01, fil_path='./', scale=1e-19, f_chind=True, figpdf=Fal
     plots
 
     '''
-    # if f_silence:
-    #     try:
-    #         import matplotlib
-    #         matplotlib.use("Agg")
-    #     except:
-    #         print('matplotlib Agg backend is not found.')
-    #         pass
-    # else:
-    #     try:
-    #         matplotlib.use("MacOSX")
-    #     except:
-    #         pass
-
-    def gaus(x,a,x0,sigma):
-        return a*exp(-(x-x0)**2/(2*sigma**2))
-
-    print('\n### Running plot_sed ###\n')
-
-    lcb = '#4682b4' # line color, blue
+    MB.logger.info('Running plot_sed')
 
     fnc  = MB.fnc 
     bfnc = MB.bfnc
@@ -132,7 +114,7 @@ def gaus(x,a,x0,sigma):
     M50 = hdul[1].data['ms'][1]
     M84 = hdul[1].data['ms'][2]
     if verbose:
-        print('Total stellar mass is %.2e'%(M50))
+        MB.logger.info('Total stellar mass is %.2e'%(M50))
 
     # Amplitude MC
     A50 = np.zeros(len(age), dtype='float')
@@ -186,7 +168,7 @@ def gaus(x,a,x0,sigma):
         DFILT = [x.strip() for x in DFILT.split(',')]
         DFWFILT = fil_fwhm(DFILT, DIR_FILT)
         if verbose:
-            print('Total dust mass is %.2e'%(MD50))
+            MB.logger.info('Total dust mass is %.2e'%(MD50))
 
     chi = hdul[1].data['chi'][0]
     chin = hdul[1].data['chi'][1]
@@ -286,7 +268,7 @@ def gaus(x,a,x0,sigma):
             ax3t = ax1.inset_axes((0.7,.35,.28,.25))
 
         f_plot_resid = False
-        print('Grism data. f_plot_resid is turned off.')
+        MB.logger.info('Grism data. f_plot_resid is turned off.')
     else:
         if f_plot_resid:
             fig_mosaic = """
@@ -312,8 +294,9 @@ def gaus(x,a,x0,sigma):
             scale = 10**(int(np.log10(np.nanmax(fybb[conbb_hs] * c / np.square(xbb[conbb_hs])) / MB.d))) / 10
         else:
             scale = 1e-19
-            print('no data point has SN > %.1f. Setting scale to %.1e'%(SNlim, scale))
-    d = MB.d * scale
+            MB.logger.info('no data point has SN > %.1f. Setting scale to %.1e'%(SNlim, scale))
+    d_scale = MB.d * scale
+    d = d_scale
 
     #######################################
     # D.Kelson like Box for BB photometry
@@ -323,23 +306,23 @@ def gaus(x,a,x0,sigma):
         for ii in range(len(xbb)):
             if eybb[ii]<100 and fybb[ii]/eybb[ii]>1:
                 xx = [xbb[ii]-exbb[ii],xbb[ii]-exbb[ii]]
-                yy = [(fybb[ii]-eybb[ii])*c/np.square(xbb[ii])/d, (fybb[ii]+eybb[ii])*c/np.square(xbb[ii])/d]
+                yy = [(fybb[ii]-eybb[ii])*c/np.square(xbb[ii])/d_scale, (fybb[ii]+eybb[ii])*c/np.square(xbb[ii])/d_scale]
                 ax1.plot(xx, yy, color='k', linestyle='-', linewidth=0.5, zorder=3)
                 xx = [xbb[ii]+exbb[ii],xbb[ii]+exbb[ii]]
-                yy = [(fybb[ii]-eybb[ii])*c/np.square(xbb[ii])/d, (fybb[ii]+eybb[ii])*c/np.square(xbb[ii])/d]
+                yy = [(fybb[ii]-eybb[ii])*c/np.square(xbb[ii])/d_scale, (fybb[ii]+eybb[ii])*c/np.square(xbb[ii])/d_scale]
                 ax1.plot(xx, yy, color='k', linestyle='-', linewidth=0.5, zorder=3)
                 xx = [xbb[ii]-exbb[ii],xbb[ii]+exbb[ii]]
-                yy = [(fybb[ii]-eybb[ii])*c/np.square(xbb[ii])/d, (fybb[ii]-eybb[ii])*c/np.square(xbb[ii])/d]
+                yy = [(fybb[ii]-eybb[ii])*c/np.square(xbb[ii])/d_scale, (fybb[ii]-eybb[ii])*c/np.square(xbb[ii])/d_scale]
                 ax1.plot(xx, yy, color='k', linestyle='-', linewidth=0.5, zorder=3)
                 xx = [xbb[ii]-exbb[ii],xbb[ii]+exbb[ii]]
-                yy = [(fybb[ii]+eybb[ii])*c/np.square(xbb[ii])/d, (fybb[ii]+eybb[ii])*c/np.square(xbb[ii])/d]
+                yy = [(fybb[ii]+eybb[ii])*c/np.square(xbb[ii])/d_scale, (fybb[ii]+eybb[ii])*c/np.square(xbb[ii])/d_scale]
                 ax1.plot(xx, yy, color='k', linestyle='-', linewidth=0.5, zorder=3)
     else: # Normal BB plot;
         # Detection;
         conbb_hs = (fybb/eybb>SNlim)
-        ax1.errorbar(xbb[conbb_hs], fybb[conbb_hs] * c / np.square(xbb[conbb_hs]) / d, \
-        yerr=eybb[conbb_hs]*c/np.square(xbb[conbb_hs])/d, color='k', linestyle='', linewidth=0.5, zorder=4)
-        ax1.plot(xbb[conbb_hs], fybb[conbb_hs] * c / np.square(xbb[conbb_hs]) / d, \
+        ax1.errorbar(xbb[conbb_hs], fybb[conbb_hs] * c / np.square(xbb[conbb_hs]) /d_scale, \
+        yerr=eybb[conbb_hs]*c/np.square(xbb[conbb_hs])/d_scale, color='k', linestyle='', linewidth=0.5, zorder=4)
+        ax1.plot(xbb[conbb_hs], fybb[conbb_hs] * c / np.square(xbb[conbb_hs]) /d_scale, \
         marker='.', color=col_dat, linestyle='', linewidth=0, zorder=4, ms=8)#, label='Obs.(BB)')
         try:
             # For any data removed fron fit (i.e. IRAC excess):
@@ -351,7 +334,7 @@ def gaus(x,a,x0,sigma):
         # Upperlim;
         sigma = 1.0
         if len(fybb[conbb_hs]):
-            leng = np.nanmax(fybb[conbb_hs] * c / np.square(xbb[conbb_hs]) / d) * 0.05 #0.2
+            leng = np.nanmax(fybb[conbb_hs] * c / np.square(xbb[conbb_hs]) /d_scale) * 0.05 #0.2
         else:
             leng = None
         conebb_ls = (fybb/eybb<=SNlim) & (eybb>0)
@@ -360,8 +343,8 @@ def gaus(x,a,x0,sigma):
             if NRbb[ii] in NR_ex[:]:
                 conebb_ls[ii] = False
 
-        ax1.errorbar(xbb[conebb_ls], eybb[conebb_ls] * c / np.square(xbb[conebb_ls]) / d * sigma, yerr=leng,\
-            uplims=eybb[conebb_ls] * c / np.square(xbb[conebb_ls]) / d * sigma, linestyle='', color=col_dat, marker='', ms=4, label='', zorder=4, capsize=3)
+        ax1.errorbar(xbb[conebb_ls], eybb[conebb_ls] * c / np.square(xbb[conebb_ls]) /d_scale * sigma, yerr=leng,\
+            uplims=eybb[conebb_ls] * c / np.square(xbb[conebb_ls]) /d_scale * sigma, linestyle='', color=col_dat, marker='', ms=4, label='', zorder=4, capsize=3)
 
 
     # For any data removed fron fit (i.e. IRAC excess):
@@ -377,11 +360,11 @@ def gaus(x,a,x0,sigma):
         x_ex, fy_ex, ey_ex, ex_ex = MB.data['bb_obs_removed']['x'], MB.data['bb_obs_removed']['fy'], MB.data['bb_obs_removed']['ey'], MB.data['bb_obs_removed']['ex']
 
         ax1.errorbar(
-            x_ex, fy_ex * c / np.square(x_ex) / d,
-            xerr=ex_ex, yerr=ey_ex*c/np.square(x_ex)/d, color='k', linestyle='', linewidth=0.5, zorder=5
+            x_ex, fy_ex * c / np.square(x_ex) /d_scale,
+            xerr=ex_ex, yerr=ey_ex*c/np.square(x_ex)/d_scale, color='k', linestyle='', linewidth=0.5, zorder=5
             )
         ax1.scatter(
-            x_ex, fy_ex * c / np.square(x_ex) / d, marker='s', color=col_ex, edgecolor='k', zorder=5, s=30
+            x_ex, fy_ex * c / np.square(x_ex) /d_scale, marker='s', color=col_ex, edgecolor='k', zorder=5, s=30
             )
         f_exclude = True
     except:
@@ -424,23 +407,23 @@ def gaus(x,a,x0,sigma):
 
         try:
             conbbd_hs = (fybbd/eybbd>SNlim)
-            ax1.errorbar(xbbd[conbbd_hs], fybbd[conbbd_hs] * c / np.square(xbbd[conbbd_hs]) / d, \
-            yerr=eybbd[conbbd_hs]*c/np.square(xbbd[conbbd_hs])/d, color='k', linestyle='', linewidth=0.5, zorder=4)
-            ax1.plot(xbbd[conbbd_hs], fybbd[conbbd_hs] * c / np.square(xbbd[conbbd_hs]) / d, \
+            ax1.errorbar(xbbd[conbbd_hs], fybbd[conbbd_hs] * c / np.square(xbbd[conbbd_hs]) /d_scale, \
+            yerr=eybbd[conbbd_hs]*c/np.square(xbbd[conbbd_hs])/d_scale, color='k', linestyle='', linewidth=0.5, zorder=4)
+            ax1.plot(xbbd[conbbd_hs], fybbd[conbbd_hs] * c / np.square(xbbd[conbbd_hs]) /d_scale, \
             '.r', linestyle='', linewidth=0, zorder=4)#, label='Obs.(BB)')
-            ax3t.plot(xbbd[conbbd_hs], fybbd[conbbd_hs] * c / np.square(xbbd[conbbd_hs]) / d, \
+            ax3t.plot(xbbd[conbbd_hs], fybbd[conbbd_hs] * c / np.square(xbbd[conbbd_hs]) /d_scale, \
             '.r', linestyle='', linewidth=0, zorder=4)#, label='Obs.(BB)')
         except:
             pass
 
         try:
             conebbd_ls = (fybbd/eybbd<=SNlim)
-            ax1.errorbar(xbbd[conebbd_ls], eybbd[conebbd_ls] * c / np.square(xbbd[conebbd_ls]) / d, \
-            yerr=fybbd[conebbd_ls]*0+np.max(fybbd[conebbd_ls]*c/np.square(xbbd[conebbd_ls])/d)*0.05, \
-            uplims=eybbd[conebbd_ls]*c/np.square(xbbd[conebbd_ls])/d, color='r', linestyle='', linewidth=0.5, zorder=4)
-            ax3t.errorbar(xbbd[conebbd_ls], eybbd[conebbd_ls] * c / np.square(xbbd[conebbd_ls]) / d, \
-            yerr=fybbd[conebbd_ls]*0+np.max(fybbd[conebbd_ls]*c/np.square(xbbd[conebbd_ls])/d)*0.05, \
-            uplims=eybbd[conebbd_ls]*c/np.square(xbbd[conebbd_ls])/d, color='r', linestyle='', linewidth=0.5, zorder=4)
+            ax1.errorbar(xbbd[conebbd_ls], eybbd[conebbd_ls] * c / np.square(xbbd[conebbd_ls]) /d_scale, \
+            yerr=fybbd[conebbd_ls]*0+np.max(fybbd[conebbd_ls]*c/np.square(xbbd[conebbd_ls])/d_scale)*0.05, \
+            uplims=eybbd[conebbd_ls]*c/np.square(xbbd[conebbd_ls])/d_scale, color='r', linestyle='', linewidth=0.5, zorder=4)
+            ax3t.errorbar(xbbd[conebbd_ls], eybbd[conebbd_ls] * c / np.square(xbbd[conebbd_ls]) /d_scale, \
+            yerr=fybbd[conebbd_ls]*0+np.max(fybbd[conebbd_ls]*c/np.square(xbbd[conebbd_ls])/d_scale)*0.05, \
+            uplims=eybbd[conebbd_ls]*c/np.square(xbbd[conebbd_ls])/d_scale, color='r', linestyle='', linewidth=0.5, zorder=4)
         except:
             pass
 
@@ -461,13 +444,13 @@ def gaus(x,a,x0,sigma):
 
             f_50_comp = np.zeros((len(age),len(y0)),'float') 
             # Keep each component;
-            f_50_comp[ii,:] = y0[:] * c / np.square(x0) / d
+            f_50_comp[ii,:] = y0[:] * c / np.square(x0) /d_scale
 
             if MB.f_dust:
                 ysump[:] += y0d_cut[:nopt]
                 ysump = np.append(ysump,y0d_cut[nopt:])
                 # Keep each component;
-                f_50_comp_dust = y0d * c / np.square(x0d) / d
+                f_50_comp_dust = y0d * c / np.square(x0d) /d_scale
 
             if MB.fneb: 
                 # Only at one age pixel;
@@ -482,7 +465,7 @@ def gaus(x,a,x0,sigma):
             ysum += y0_r
             ysump[:nopt] += y0p
 
-            f_50_comp[ii,:] = y0_r[:] * c / np.square(x0_tmp) / d
+            f_50_comp[ii,:] = y0_r[:] * c / np.square(x0_tmp) /d_scale
 
         # The following needs revised.
         f_uvj = False
@@ -516,9 +499,9 @@ def gaus(x,a,x0,sigma):
     #############
     conbb_ymax = (xbb>0) & (fybb>0) & (eybb>0) & (fybb/eybb>SNlim)
     if len(fybb[conbb_ymax]):
-        ymax = np.nanmax(fybb[conbb_ymax]*c/np.square(xbb[conbb_ymax])/d) * 1.6
+        ymax = np.nanmax(fybb[conbb_ymax]*c/np.square(xbb[conbb_ymax])/d_scale) * 1.6
     else:
-        ymax = np.nanmax(fybb*c/np.square(xbb)/d) * 1.6
+        ymax = np.nanmax(fybb*c/np.square(xbb)/d_scale) * 1.6
 
     ax1.set_xlabel('Observed wavelength [$\mathrm{\mu m}$]', fontsize=11)
     ax1.set_ylabel('$f_\lambda$ [$10^{%d}\mathrm{erg}/\mathrm{s}/\mathrm{cm}^{2}/\mathrm{\AA}$]'%(np.log10(scale)),fontsize=11,labelpad=2)
@@ -612,7 +595,7 @@ def gaus(x,a,x0,sigma):
     DL = MB.cosmo.luminosity_distance(zbes).value * Mpc_cm #, **cosmo) # Luminositydistance in cm
     Cons = (4.*np.pi*DL**2/(1.+zbes))
     if f_grsm:
-        print('This function (write_lines) needs to be revised.')
+        MB.logger.warning('This function (write_lines) needs to be revised.')
         write_lines(ID, zbes, DIR_OUT=MB.DIR_OUT)
 
 
@@ -620,9 +603,9 @@ def gaus(x,a,x0,sigma):
     # Zoom in Line regions
     ##########################
     if f_grsm:
-        ax2t.errorbar(xg2, fg2 * c/np.square(xg2)/d, yerr=eg2 * c/np.square(xg2)/d, lw=0.5, color='#DF4E00', zorder=10, alpha=1., label='', capsize=0)
-        ax2t.errorbar(xg1, fg1 * c/np.square(xg1)/d, yerr=eg1 * c/np.square(xg1)/d, lw=0.5, color='g', zorder=10, alpha=1., label='', capsize=0)
-        ax2t.errorbar(xg0, fg0 * c/np.square(xg0)/d, yerr=eg0 * c/np.square(xg0)/d, lw=0.5, color='royalblue', zorder=10, alpha=1., label='', capsize=0)
+        ax2t.errorbar(xg2, fg2 * c/np.square(xg2)/d_scale, yerr=eg2 * c/np.square(xg2)/d_scale, lw=0.5, color='#DF4E00', zorder=10, alpha=1., label='', capsize=0)
+        ax2t.errorbar(xg1, fg1 * c/np.square(xg1)/d_scale, yerr=eg1 * c/np.square(xg1)/d_scale, lw=0.5, color='g', zorder=10, alpha=1., label='', capsize=0)
+        ax2t.errorbar(xg0, fg0 * c/np.square(xg0)/d_scale, yerr=eg0 * c/np.square(xg0)/d_scale, lw=0.5, color='royalblue', zorder=10, alpha=1., label='', capsize=0)
 
         xgrism = np.concatenate([xg0,xg1,xg2])
         fgrism = np.concatenate([fg0,fg1,fg2])
@@ -727,7 +710,7 @@ def gaus(x,a,x0,sigma):
                 fm_tmp_nl += mod0_tmp
 
             # Each;
-            ytmp_each[kk,:,ss] = mod0_tmp[:] * c / np.square(xm_tmp[:]) / d
+            ytmp_each[kk,:,ss] = mod0_tmp[:] * c / np.square(xm_tmp[:]) /d_scale
 
         #
         # Dust component;
@@ -758,25 +741,25 @@ def gaus(x,a,x0,sigma):
                 ytmp_dust = np.zeros((mmax,len(x1_dust)), dtype='float')
                 ytmp_comp = np.zeros((mmax,len(x1_tot)), dtype='float')
 
-            ytmp_dust[kk,:] = model_dust * c/np.square(x1_dust)/d
+            ytmp_dust[kk,:] = model_dust * c/np.square(x1_dust)/d_scale
             model_tot = np.interp(x1_tot,xx_tmp,fm_tmp) + np.interp(x1_tot,x1_dust,model_dust)
             model_tot_nl = np.interp(x1_tot,xx_tmp,fm_tmp_nl) + np.interp(x1_tot,x1_dust,model_dust)
 
-            ytmp[kk,:] = model_tot[:] * c/np.square(x1_tot[:])/d
-            ytmp_nl[kk,:] = model_tot_nl[:] * c/np.square(x1_tot[:])/d
+            ytmp[kk,:] = model_tot[:] * c/np.square(x1_tot[:])/d_scale
+            ytmp_nl[kk,:] = model_tot_nl[:] * c/np.square(x1_tot[:])/d_scale
 
         else:
             x1_tot = xm_tmp
-            ytmp[kk,:] = fm_tmp[:] * c / np.square(xm_tmp[:]) / d
-            ytmp_nl[kk,:] = fm_tmp_nl[:] * c / np.square(xm_tmp[:]) / d
+            ytmp[kk,:] = fm_tmp[:] * c / np.square(xm_tmp[:]) /d_scale
+            ytmp_nl[kk,:] = fm_tmp_nl[:] * c / np.square(xm_tmp[:]) /d_scale
 
         # Get FUV flux density;
-        Fuv[kk] = get_Fuv(x1_tot[:]/(1.+zbes), (ytmp[kk,:]/(c/np.square(x1_tot)/d)) * (DL**2/(1.+zbes)) / (DL10**2), lmin=1250, lmax=1650)
-        Fuv28[kk] = get_Fuv(x1_tot[:]/(1.+zbes), (ytmp[kk,:]/(c/np.square(x1_tot)/d)) * (4*np.pi*DL**2/(1.+zbes))*Cmznu, lmin=1500, lmax=2800)
+        Fuv[kk] = get_Fuv(x1_tot[:]/(1.+zbes), (ytmp[kk,:]/(c/np.square(x1_tot)/d_scale)) * (DL**2/(1.+zbes)) / (DL10**2), lmin=1250, lmax=1650)
+        Fuv28[kk] = get_Fuv(x1_tot[:]/(1.+zbes), (ytmp[kk,:]/(c/np.square(x1_tot)/d_scale)) * (4*np.pi*DL**2/(1.+zbes))*Cmznu, lmin=1500, lmax=2800)
         Lir[kk] = 0
 
         # Get UVJ Color;
-        _,fconv = filconv_fast(MB.filts_rf, MB.band_rf, x1_tot[:]/(1.+zbes), (ytmp[kk,:]/(c/np.square(x1_tot)/d)))
+        _,fconv = filconv_fast(MB.filts_rf, MB.band_rf, x1_tot[:]/(1.+zbes), (ytmp[kk,:]/(c/np.square(x1_tot)/d_scale)))
         UVJ[kk,0] = -2.5*np.log10(fconv[0]/fconv[2])
         UVJ[kk,1] = -2.5*np.log10(fconv[1]/fconv[2])
         UVJ[kk,2] = -2.5*np.log10(fconv[2]/fconv[3])
@@ -823,8 +806,8 @@ def gaus(x,a,x0,sigma):
         
     # Attach the data point in MB;
     MB.sed_wave_obs = xbb
-    MB.sed_flux_obs = fybb * c / np.square(xbb) / d
-    MB.sed_eflux_obs = eybb * c / np.square(xbb) / d
+    MB.sed_flux_obs = fybb * c / np.square(xbb) /d_scale
+    MB.sed_eflux_obs = eybb * c / np.square(xbb) /d_scale
     # Attach the best SED to MB;
     MB.sed_wave = x1_tot
     MB.sed_flux16 = ytmp16
@@ -847,7 +830,7 @@ def gaus(x,a,x0,sigma):
                 ax1.plot(x1_tot[::nstep_plot], np.percentile(ysumtmp2[:, ::nstep_plot], 50, axis=0), linestyle='--', lw=.5, color=col[ii], alpha=alp_fancy, zorder=-3)
             ysumtmp2_prior[:] = np.percentile(ysumtmp2[:, :], 50, axis=0)
     elif f_fill:
-        print('f_fancyplot is False. f_fill is set to False.')
+        MB.logger.info('f_fancyplot is False. f_fill is set to False.')
 
     # Calculate non-det chi2
     # based on Sawick12
@@ -930,10 +913,10 @@ def gaus(x,a,x0,sigma):
 
         if f_plot_resid:
             conbb_hs = (fybb/eybb>SNlim)
-            axes['B'].scatter(lbb[iix][conbb_hs], ((fybb*c/np.square(xbb)/d - fbb)/(eybb*c/np.square(xbb)/d))[iix][conbb_hs], lw=0.5, color='none', edgecolor='r', zorder=3, alpha=1.0, marker='.', s=50)
+            axes['B'].scatter(lbb[iix][conbb_hs], ((fybb*c/np.square(xbb)/d_scale - fbb)/(eybb*c/np.square(xbb)/d_scale))[iix][conbb_hs], lw=0.5, color='none', edgecolor='r', zorder=3, alpha=1.0, marker='.', s=50)
             conbb_hs = (fybb/eybb<=SNlim) & (eybb>0)
-            axes['B'].errorbar(lbb[iix][conbb_hs], ((eybb*c/np.square(xbb)/d - fbb)/(eybb*c/np.square(xbb)/d))[iix][conbb_hs], yerr=leng,\
-                uplims=((fybb*c/np.square(xbb)/d - fbb)/(eybb*c/np.square(xbb)/d))[iix][conbb_hs] * sigma, linestyle='',\
+            axes['B'].errorbar(lbb[iix][conbb_hs], ((eybb*c/np.square(xbb)/d_scale - fbb)/(eybb*c/np.square(xbb)/d_scale))[iix][conbb_hs], yerr=leng,\
+                uplims=((fybb*c/np.square(xbb)/d_scale - fbb)/(eybb*c/np.square(xbb)/d_scale))[iix][conbb_hs] * sigma, linestyle='',\
                 color=col_dat, lw=0.5, marker='', ms=4, label='', zorder=4, capsize=1.5)
             axes['B'].set_xscale(ax1.get_xscale())
             axes['B'].set_xlim(ax1.get_xlim())
@@ -955,12 +938,12 @@ def gaus(x,a,x0,sigma):
 
             # Rest-frame EW;
             # Note about 16/84 in fbb
-            EW16 = (fy_ex * c / np.square(x_ex) / d - fbb84[iix2]) / (fbb[iix2]) * lfwhm[iix2] / (1.+zbes)
-            EW50 = (fy_ex * c / np.square(x_ex) / d - fbb[iix2]) / (fbb[iix2]) * lfwhm[iix2] / (1.+zbes)
-            EW84 = (fy_ex * c / np.square(x_ex) / d - fbb16[iix2]) / (fbb[iix2]) * lfwhm[iix2] / (1.+zbes)
+            EW16 = (fy_ex * c / np.square(x_ex) /d_scale - fbb84[iix2]) / (fbb[iix2]) * lfwhm[iix2] / (1.+zbes)
+            EW50 = (fy_ex * c / np.square(x_ex) /d_scale - fbb[iix2]) / (fbb[iix2]) * lfwhm[iix2] / (1.+zbes)
+            EW84 = (fy_ex * c / np.square(x_ex) /d_scale - fbb16[iix2]) / (fbb[iix2]) * lfwhm[iix2] / (1.+zbes)
 
-            EW50_er1 = ((fy_ex-ey_ex) * c / np.square(x_ex) / d - fbb[iix2]) / (fbb[iix2]) * lfwhm[iix2] / (1.+zbes)
-            EW50_er2 = ((fy_ex+ey_ex) * c / np.square(x_ex) / d - fbb[iix2]) / (fbb[iix2]) * lfwhm[iix2] / (1.+zbes)
+            EW50_er1 = ((fy_ex-ey_ex) * c / np.square(x_ex) /d_scale - fbb[iix2]) / (fbb[iix2]) * lfwhm[iix2] / (1.+zbes)
+            EW50_er2 = ((fy_ex+ey_ex) * c / np.square(x_ex) /d_scale - fbb[iix2]) / (fbb[iix2]) * lfwhm[iix2] / (1.+zbes)
 
             cnt50 = fbb[iix2]
             cnt16 = fbb16[iix2]
@@ -1040,9 +1023,9 @@ def gaus(x,a,x0,sigma):
 
         col5 = fits.Column(name='wave_obs', format='E', unit='AA', array=xbb)
         col00.append(col5)
-        col6 = fits.Column(name='f_obs', format='E', unit='1e%derg/s/cm2/AA'%(np.log10(scale)), array=fybb[:] * c / np.square(xbb[:]) / d)
+        col6 = fits.Column(name='f_obs', format='E', unit='1e%derg/s/cm2/AA'%(np.log10(scale)), array=fybb[:] * c / np.square(xbb[:]) /d_scale)
         col00.append(col6)
-        col7 = fits.Column(name='e_obs', format='E', unit='1e%derg/s/cm2/AA'%(np.log10(scale)), array=eybb[:] * c / np.square(xbb[:]) / d)
+        col7 = fits.Column(name='e_obs', format='E', unit='1e%derg/s/cm2/AA'%(np.log10(scale)), array=eybb[:] * c / np.square(xbb[:]) /d_scale)
         col00.append(col7)
 
         hdr = fits.Header()
@@ -1186,9 +1169,9 @@ def gaus(x,a,x0,sigma):
                 tree_spec['model'].update({'EW_%s_84'%(ew_label[ii]): EW84[ii] * u.AA})
                 tree_spec['model'].update({'EW_%s_e1'%(ew_label[ii]): EW50_er1[ii] * u.AA})
                 tree_spec['model'].update({'EW_%s_e2'%(ew_label[ii]): EW50_er2[ii] * u.AA})
-                tree_spec['model'].update({'cnt_%s_16'%(ew_label[ii]): cnt16[ii] / (c / np.square(x1_tot[:]) / d) * Cnu_to_Jy * u.uJy})
-                tree_spec['model'].update({'cnt_%s_50'%(ew_label[ii]): cnt50[ii] / (c / np.square(x1_tot[:]) / d) * Cnu_to_Jy * u.uJy})
-                tree_spec['model'].update({'cnt_%s_84'%(ew_label[ii]): cnt84[ii] / (c / np.square(x1_tot[:]) / d) * Cnu_to_Jy * u.uJy})
+                tree_spec['model'].update({'cnt_%s_16'%(ew_label[ii]): cnt16[ii] / (c / np.square(x1_tot[:]) /d_scale) * Cnu_to_Jy * u.uJy})
+                tree_spec['model'].update({'cnt_%s_50'%(ew_label[ii]): cnt50[ii] / (c / np.square(x1_tot[:]) /d_scale) * Cnu_to_Jy * u.uJy})
+                tree_spec['model'].update({'cnt_%s_84'%(ew_label[ii]): cnt84[ii] / (c / np.square(x1_tot[:]) /d_scale) * Cnu_to_Jy * u.uJy})
                 tree_spec['model'].update({'L_%s_16'%(ew_label[ii]): L16[ii] * u.erg / u.s})
                 tree_spec['model'].update({'L_%s_50'%(ew_label[ii]): L50[ii] * u.erg / u.s})
                 tree_spec['model'].update({'L_%s_84'%(ew_label[ii]): L84[ii] * u.erg / u.s})
@@ -1261,8 +1244,8 @@ def gaus(x,a,x0,sigma):
             ax2t.set_xlim(xgmin, xgmax)
 
             conaa = (x0>xgmin-50) & (x0<xgmax+50)
-            ymaxzoom = np.max(ysum[conaa]*c/np.square(x0[conaa])/d) * 1.15
-            yminzoom = np.min(ysum[conaa]*c/np.square(x0[conaa])/d) / 1.15
+            ymaxzoom = np.max(ysum[conaa]*c/np.square(x0[conaa])/d_scale) * 1.15
+            yminzoom = np.min(ysum[conaa]*c/np.square(x0[conaa])/d_scale) / 1.15
 
             ax2t.set_ylim(yminzoom, ymaxzoom)
             ax2t.xaxis.labelpad = -2
@@ -1280,8 +1263,8 @@ def gaus(x,a,x0,sigma):
             ax2t.set_xlim(xgmin, xgmax)
 
             conaa = (x0>xgmin-50) & (x0<xgmax+50)
-            ymaxzoom = np.max(ysum[conaa]*c/np.square(x0[conaa])/d) * 1.15
-            yminzoom = np.min(ysum[conaa]*c/np.square(x0[conaa])/d) / 1.15
+            ymaxzoom = np.max(ysum[conaa]*c/np.square(x0[conaa])/d_scale) * 1.15
+            yminzoom = np.min(ysum[conaa]*c/np.square(x0[conaa])/d_scale) / 1.15
 
             ax2t.set_ylim(yminzoom, ymaxzoom)
             ax2t.xaxis.labelpad = -2
@@ -1295,7 +1278,7 @@ def gaus(x,a,x0,sigma):
     if MB.f_dust:
         try:
             contmp = (x1_tot>10*1e4) #& (fybbd/eybbd>SNlim)
-            y3min, y3max = -.2*np.max((model_tot * c/ np.square(x1_tot) / d)[contmp]), np.max((model_tot * c/ np.square(x1_tot) / d)[contmp])*2.0
+            y3min, y3max = -.2*np.max((model_tot * c/ np.square(x1_tot) /d_scale)[contmp]), np.max((model_tot * c/ np.square(x1_tot) /d_scale)[contmp])*2.0
             ax3t.set_ylim(y3min, y3max)
         except:
             if verbose:
@@ -1666,23 +1649,23 @@ def gaus(x,a,x0,sigma):
         for ii in range(len(xbb)):
             if eybb[ii]<100 and fybb[ii]/eybb[ii]>1:
                 xx = [xbb[ii]-exbb[ii],xbb[ii]-exbb[ii]]
-                yy = [(fybb[ii]-eybb[ii])*c/np.square(xbb[ii])/d, (fybb[ii]+eybb[ii])*c/np.square(xbb[ii])/d]
+                yy = [(fybb[ii]-eybb[ii])*c/np.square(xbb[ii])/d_scale, (fybb[ii]+eybb[ii])*c/np.square(xbb[ii])/d_scale]
                 ax1.plot(xx, yy, color='k', linestyle='-', linewidth=0.5, zorder=3)
                 xx = [xbb[ii]+exbb[ii],xbb[ii]+exbb[ii]]
-                yy = [(fybb[ii]-eybb[ii])*c/np.square(xbb[ii])/d, (fybb[ii]+eybb[ii])*c/np.square(xbb[ii])/d]
+                yy = [(fybb[ii]-eybb[ii])*c/np.square(xbb[ii])/d_scale, (fybb[ii]+eybb[ii])*c/np.square(xbb[ii])/d_scale]
                 ax1.plot(xx, yy, color='k', linestyle='-', linewidth=0.5, zorder=3)
                 xx = [xbb[ii]-exbb[ii],xbb[ii]+exbb[ii]]
-                yy = [(fybb[ii]-eybb[ii])*c/np.square(xbb[ii])/d, (fybb[ii]-eybb[ii])*c/np.square(xbb[ii])/d]
+                yy = [(fybb[ii]-eybb[ii])*c/np.square(xbb[ii])/d_scale, (fybb[ii]-eybb[ii])*c/np.square(xbb[ii])/d_scale]
                 ax1.plot(xx, yy, color='k', linestyle='-', linewidth=0.5, zorder=3)
                 xx = [xbb[ii]-exbb[ii],xbb[ii]+exbb[ii]]
-                yy = [(fybb[ii]+eybb[ii])*c/np.square(xbb[ii])/d, (fybb[ii]+eybb[ii])*c/np.square(xbb[ii])/d]
+                yy = [(fybb[ii]+eybb[ii])*c/np.square(xbb[ii])/d_scale, (fybb[ii]+eybb[ii])*c/np.square(xbb[ii])/d_scale]
                 ax1.plot(xx, yy, color='k', linestyle='-', linewidth=0.5, zorder=3)
     else: # Normal BB plot;
         # Detection;
         conbb_hs = (fybb/eybb>SNlim)
-        ax1.errorbar(xbb[conbb_hs], fybb[conbb_hs] * c / np.square(xbb[conbb_hs]) / d, \
-        yerr=eybb[conbb_hs]*c/np.square(xbb[conbb_hs])/d, color='k', linestyle='', linewidth=0.5, zorder=4)
-        ax1.plot(xbb[conbb_hs], fybb[conbb_hs] * c / np.square(xbb[conbb_hs]) / d, \
+        ax1.errorbar(xbb[conbb_hs], fybb[conbb_hs] * c / np.square(xbb[conbb_hs]) /d_scale, \
+        yerr=eybb[conbb_hs]*c/np.square(xbb[conbb_hs])/d_scale, color='k', linestyle='', linewidth=0.5, zorder=4)
+        ax1.plot(xbb[conbb_hs], fybb[conbb_hs] * c / np.square(xbb[conbb_hs]) /d_scale, \
         marker='.', color=col_dat, linestyle='', linewidth=0, zorder=4, ms=8)#, label='Obs.(BB)')
 
         try:
@@ -1694,15 +1677,15 @@ def gaus(x,a,x0,sigma):
 
         # Upperlim;
         sigma = 1.0
-        leng = np.max(fybb[conbb_hs] * c / np.square(xbb[conbb_hs]) / d) * 0.05 #0.2
+        leng = np.max(fybb[conbb_hs] * c / np.square(xbb[conbb_hs]) /d_scale) * 0.05 #0.2
         conebb_ls = (fybb/eybb<=SNlim) & (eybb>0)
         
         for ii in range(len(xbb)):
             if NR[ii] in NR_ex[:]:
                 conebb_ls[ii] = False
         
-        ax1.errorbar(xbb[conebb_ls], eybb[conebb_ls] * c / np.square(xbb[conebb_ls]) / d * sigma, yerr=leng,\
-            uplims=eybb[conebb_ls] * c / np.square(xbb[conebb_ls]) / d * sigma, linestyle='',color=col_dat, marker='', ms=4, label='', zorder=4, capsize=3)
+        ax1.errorbar(xbb[conebb_ls], eybb[conebb_ls] * c / np.square(xbb[conebb_ls]) /d_scale * sigma, yerr=leng,\
+            uplims=eybb[conebb_ls] * c / np.square(xbb[conebb_ls]) /d_scale * sigma, linestyle='',color=col_dat, marker='', ms=4, label='', zorder=4, capsize=3)
 
 
     # For any data removed fron fit (i.e. IRAC excess):
@@ -1718,9 +1701,9 @@ def gaus(x,a,x0,sigma):
         ey_ex = data_ex['col4']
         ex_ex = data_ex['col5']
 
-        ax1.errorbar(x_ex, fy_ex * c / np.square(x_ex) / d, \
-        xerr=ex_ex, yerr=ey_ex*c/np.square(x_ex)/d, color='k', linestyle='', linewidth=0.5, zorder=5)
-        ax1.scatter(x_ex, fy_ex * c / np.square(x_ex) / d, marker='s', color=col_ex, edgecolor='k', zorder=5, s=30)
+        ax1.errorbar(x_ex, fy_ex * c / np.square(x_ex) /d_scale, \
+        xerr=ex_ex, yerr=ey_ex*c/np.square(x_ex)/d_scale, color='k', linestyle='', linewidth=0.5, zorder=5)
+        ax1.scatter(x_ex, fy_ex * c / np.square(x_ex) /d_scale, marker='s', color=col_ex, edgecolor='k', zorder=5, s=30)
         f_exclude = True
     except:
         pass
@@ -1764,23 +1747,23 @@ def gaus(x,a,x0,sigma):
 
         try:
             conbbd_hs = (fybbd/eybbd>SNlim)
-            ax1.errorbar(xbbd[conbbd_hs], fybbd[conbbd_hs] * c / np.square(xbbd[conbbd_hs]) / d, \
-            yerr=eybbd[conbbd_hs]*c/np.square(xbbd[conbbd_hs])/d, color='k', linestyle='', linewidth=0.5, zorder=4)
-            ax1.plot(xbbd[conbbd_hs], fybbd[conbbd_hs] * c / np.square(xbbd[conbbd_hs]) / d, \
+            ax1.errorbar(xbbd[conbbd_hs], fybbd[conbbd_hs] * c / np.square(xbbd[conbbd_hs]) /d_scale, \
+            yerr=eybbd[conbbd_hs]*c/np.square(xbbd[conbbd_hs])/d_scale, color='k', linestyle='', linewidth=0.5, zorder=4)
+            ax1.plot(xbbd[conbbd_hs], fybbd[conbbd_hs] * c / np.square(xbbd[conbbd_hs]) /d_scale, \
             '.r', linestyle='', linewidth=0, zorder=4)#, label='Obs.(BB)')
-            ax3t.plot(xbbd[conbbd_hs], fybbd[conbbd_hs] * c / np.square(xbbd[conbbd_hs]) / d, \
+            ax3t.plot(xbbd[conbbd_hs], fybbd[conbbd_hs] * c / np.square(xbbd[conbbd_hs]) /d_scale, \
             '.r', linestyle='', linewidth=0, zorder=4)#, label='Obs.(BB)')
         except:
             pass
 
         try:
             conebbd_ls = (fybbd/eybbd<=SNlim)
-            ax1.errorbar(xbbd[conebbd_ls], eybbd[conebbd_ls] * c / np.square(xbbd[conebbd_ls]) / d, \
-            yerr=fybbd[conebbd_ls]*0+np.max(fybbd[conebbd_ls]*c/np.square(xbbd[conebbd_ls])/d)*0.05, \
-            uplims=eybbd[conebbd_ls]*c/np.square(xbbd[conebbd_ls])/d, color='r', linestyle='', linewidth=0.5, zorder=4)
-            ax3t.errorbar(xbbd[conebbd_ls], eybbd[conebbd_ls] * c / np.square(xbbd[conebbd_ls]) / d, \
-            yerr=fybbd[conebbd_ls]*0+np.max(fybbd[conebbd_ls]*c/np.square(xbbd[conebbd_ls])/d)*0.05, \
-            uplims=eybbd[conebbd_ls]*c/np.square(xbbd[conebbd_ls])/d, color='r', linestyle='', linewidth=0.5, zorder=4)
+            ax1.errorbar(xbbd[conebbd_ls], eybbd[conebbd_ls] * c / np.square(xbbd[conebbd_ls]) /d_scale, \
+            yerr=fybbd[conebbd_ls]*0+np.max(fybbd[conebbd_ls]*c/np.square(xbbd[conebbd_ls])/d_scale)*0.05, \
+            uplims=eybbd[conebbd_ls]*c/np.square(xbbd[conebbd_ls])/d_scale, color='r', linestyle='', linewidth=0.5, zorder=4)
+            ax3t.errorbar(xbbd[conebbd_ls], eybbd[conebbd_ls] * c / np.square(xbbd[conebbd_ls]) /d_scale, \
+            yerr=fybbd[conebbd_ls]*0+np.max(fybbd[conebbd_ls]*c/np.square(xbbd[conebbd_ls])/d_scale)*0.05, \
+            uplims=eybbd[conebbd_ls]*c/np.square(xbbd[conebbd_ls])/d_scale, color='r', linestyle='', linewidth=0.5, zorder=4)
         except:
             pass
 
@@ -1796,7 +1779,7 @@ def gaus(x,a,x0,sigma):
     y0, x0 = MB.fnc.get_template(vals, f_val=False, check_bound=False, lib_all=True)
 
     ysum = y0
-    f_50_comp = y0[:] * c / np.square(x0) / d
+    f_50_comp = y0[:] * c / np.square(x0) /d_scale
 
     ysump = y0p
     nopt = len(ysump)
@@ -1804,7 +1787,7 @@ def gaus(x,a,x0,sigma):
     if f_dust:
         ysump[:] += y0d_cut[:nopt]
         ysump = np.append(ysump,y0d_cut[nopt:])
-        f_50_comp_dust = y0d * c / np.square(x0d) / d
+        f_50_comp_dust = y0d * c / np.square(x0d) /d_scale
 
     if MB.fneb: 
         # Only at one age pixel;
@@ -1824,7 +1807,7 @@ def gaus(x,a,x0,sigma):
             y0keep = y0tmp
         else:
             y0keep += y0tmp
-        ax1.plot(x0tmp, y0tmp * c / np.square(x0tmp) / d, linestyle='--', lw=0.5, color=col[aa])
+        ax1.plot(x0tmp, y0tmp * c / np.square(x0tmp) /d_scale, linestyle='--', lw=0.5, color=col[aa])
         vals_each['A%d'%aa] = 0
 
     # Plot best fit;
@@ -1834,7 +1817,7 @@ def gaus(x,a,x0,sigma):
     # Main result
     #############
     conbb_ymax = (xbb>0) & (fybb>0) & (eybb>0) & (fybb/eybb>1) # (conbb) &
-    ymax = np.max(fybb[conbb_ymax]*c/np.square(xbb[conbb_ymax])/d) * 1.6
+    ymax = np.max(fybb[conbb_ymax]*c/np.square(xbb[conbb_ymax])/d_scale) * 1.6
 
     xboxl = 17000
     xboxu = 28000
@@ -1935,9 +1918,9 @@ def gaus(x,a,x0,sigma):
     ##########################
     if f_grsm:
         conspec = (NR<10000) #& (fy/ey>1)
-        ax2t.errorbar(xg2, fg2 * c/np.square(xg2)/d, yerr=eg2 * c/np.square(xg2)/d, lw=0.5, color='#DF4E00', zorder=10, alpha=1., label='', capsize=0)
-        ax2t.errorbar(xg1, fg1 * c/np.square(xg1)/d, yerr=eg1 * c/np.square(xg1)/d, lw=0.5, color='g', zorder=10, alpha=1., label='', capsize=0)
-        ax2t.errorbar(xg0, fg0 * c/np.square(xg0)/d, yerr=eg0 * c/np.square(xg0)/d, lw=0.5, linestyle='', color='royalblue', zorder=10, alpha=1., label='', capsize=0)
+        ax2t.errorbar(xg2, fg2 * c/np.square(xg2)/d_scale, yerr=eg2 * c/np.square(xg2)/d_scale, lw=0.5, color='#DF4E00', zorder=10, alpha=1., label='', capsize=0)
+        ax2t.errorbar(xg1, fg1 * c/np.square(xg1)/d_scale, yerr=eg1 * c/np.square(xg1)/d_scale, lw=0.5, color='g', zorder=10, alpha=1., label='', capsize=0)
+        ax2t.errorbar(xg0, fg0 * c/np.square(xg0)/d_scale, yerr=eg0 * c/np.square(xg0)/d_scale, lw=0.5, linestyle='', color='royalblue', zorder=10, alpha=1., label='', capsize=0)
 
         xgrism = np.concatenate([xg0,xg1,xg2])
         fgrism = np.concatenate([fg0,fg1,fg2])
@@ -2029,7 +2012,7 @@ def gaus(x,a,x0,sigma):
 
         if False:
             # Each;
-            ytmp_each[kk,:,ss] = mod0_tmp[:] * c / np.square(xm_tmp[:]) / d
+            ytmp_each[kk,:,ss] = mod0_tmp[:] * c / np.square(xm_tmp[:]) /d_scale
             #if kk == 100:
             #    ax1.plot(xm_tmp[:], ytmp_each[kk,:,ss], color=col[ss], linestyle='--')
 
@@ -2059,14 +2042,14 @@ def gaus(x,a,x0,sigma):
                 ytmp = np.zeros((mmax,len(x1_tot)), dtype='float')
                 ytmp_dust = np.zeros((mmax,len(x1_dust)), dtype='float')
 
-            ytmp_dust[kk,:] = model_dust * c/np.square(x1_dust)/d
+            ytmp_dust[kk,:] = model_dust * c/np.square(x1_dust)/d_scale
             model_tot = np.interp(x1_tot,xm_tmp,fm_tmp) + np.interp(x1_tot,x1_dust,model_dust)
 
-            ytmp[kk,:] = model_tot[:] * c/np.square(x1_tot[:])/d
+            ytmp[kk,:] = model_tot[:] * c/np.square(x1_tot[:])/d_scale
 
         else:
             x1_tot = xm_tmp
-            ytmp[kk,:] = fm_tmp[:] * c / np.square(xm_tmp[:]) / d
+            ytmp[kk,:] = fm_tmp[:] * c / np.square(xm_tmp[:]) /d_scale
 
         # plot random sed;
         plot_mc = True
@@ -2079,12 +2062,12 @@ def gaus(x,a,x0,sigma):
 
         if True:
             # Get FUV flux;
-            Fuv[kk] = get_Fuv(x1_tot[:]/(1.+zbes), (ytmp[kk,:]/(c/np.square(x1_tot)/d)) * (DL**2/(1.+zbes)) / (DL10**2), lmin=1250, lmax=1650)
-            Fuv28[kk] = get_Fuv(x1_tot[:]/(1.+zbes), (ytmp[kk,:]/(c/np.square(x1_tot)/d)) * (4*np.pi*DL**2/(1.+zbes))*Cmznu, lmin=1500, lmax=2800)
+            Fuv[kk] = get_Fuv(x1_tot[:]/(1.+zbes), (ytmp[kk,:]/(c/np.square(x1_tot)/d_scale)) * (DL**2/(1.+zbes)) / (DL10**2), lmin=1250, lmax=1650)
+            Fuv28[kk] = get_Fuv(x1_tot[:]/(1.+zbes), (ytmp[kk,:]/(c/np.square(x1_tot)/d_scale)) * (4*np.pi*DL**2/(1.+zbes))*Cmznu, lmin=1500, lmax=2800)
             Lir[kk] = 0
 
             # Get UVJ Color;
-            lmconv,fconv = filconv_fast(MB.filts_rf, MB.band_rf, x1_tot[:]/(1.+zbes), (ytmp[kk,:]/(c/np.square(x1_tot)/d)))
+            lmconv,fconv = filconv_fast(MB.filts_rf, MB.band_rf, x1_tot[:]/(1.+zbes), (ytmp[kk,:]/(c/np.square(x1_tot)/d_scale)))
             UVJ[kk,0] = -2.5*np.log10(fconv[0]/fconv[2])
             UVJ[kk,1] = -2.5*np.log10(fconv[1]/fconv[2])
             UVJ[kk,2] = -2.5*np.log10(fconv[2]/fconv[3])
@@ -2127,8 +2110,8 @@ def gaus(x,a,x0,sigma):
 
     # Attach the data point in MB;
     MB.sed_wave_obs = xbb
-    MB.sed_flux_obs = fybb * c / np.square(xbb) / d
-    MB.sed_eflux_obs = eybb * c / np.square(xbb) / d
+    MB.sed_flux_obs = fybb * c / np.square(xbb) /d_scale
+    MB.sed_eflux_obs = eybb * c / np.square(xbb) /d_scale
     # Attach the best SED to MB;
     MB.sed_wave = x1_tot
     MB.sed_flux16 = ytmp16
@@ -2220,7 +2203,7 @@ def func_tmp(xint,eobs,fmodel):
     if f_dust:
         ALLFILT = np.append(SFILT,DFILT)
         #for ii in range(len(x1_tot)):
-        #    print(x1_tot[ii], model_tot[ii]*c/np.square(x1_tot[ii])/d)
+        #    print(x1_tot[ii], model_tot[ii]*c/np.square(x1_tot[ii])/d_scale)
         lbb, fbb, lfwhm = filconv(ALLFILT, x1_tot, ytmp50, DIR_FILT, fw=True)
         lbb, fbb16, lfwhm = filconv(ALLFILT, x1_tot, ytmp16, DIR_FILT, fw=True)
         lbb, fbb84, lfwhm = filconv(ALLFILT, x1_tot, ytmp84, DIR_FILT, fw=True)
@@ -2251,10 +2234,10 @@ def func_tmp(xint,eobs,fmodel):
 
         if f_plot_resid:
             conbb_hs = (fybb/eybb>SNlim)
-            axes['B'].scatter(lbb[iix][conbb_hs], ((fybb*c/np.square(xbb)/d - fbb)/(eybb*c/np.square(xbb)/d))[iix][conbb_hs], lw=0.5, color='none', edgecolor='r', zorder=3, alpha=1.0, marker='.', s=50)
+            axes['B'].scatter(lbb[iix][conbb_hs], ((fybb*c/np.square(xbb)/d_scale - fbb)/(eybb*c/np.square(xbb)/d_scale))[iix][conbb_hs], lw=0.5, color='none', edgecolor='r', zorder=3, alpha=1.0, marker='.', s=50)
             conbb_hs = (fybb/eybb<=SNlim) & (eybb>0)
-            axes['B'].errorbar(lbb[iix][conbb_hs], ((eybb*c/np.square(xbb)/d - fbb)/(eybb*c/np.square(xbb)/d))[iix][conbb_hs], yerr=leng,\
-                uplims=((fybb*c/np.square(xbb)/d - fbb)/(eybb*c/np.square(xbb)/d))[iix][conbb_hs] * sigma, linestyle='',\
+            axes['B'].errorbar(lbb[iix][conbb_hs], ((eybb*c/np.square(xbb)/d_scale - fbb)/(eybb*c/np.square(xbb)/d_scale))[iix][conbb_hs], yerr=leng,\
+                uplims=((fybb*c/np.square(xbb)/d_scale - fbb)/(eybb*c/np.square(xbb)/d_scale))[iix][conbb_hs] * sigma, linestyle='',\
                 color=col_dat, lw=0.5, marker='', ms=4, label='', zorder=4, capsize=1.5)
             axes['B'].set_xscale(ax1.get_xscale())
             axes['B'].set_xlim(ax1.get_xlim())
@@ -2276,12 +2259,12 @@ def func_tmp(xint,eobs,fmodel):
 
             # Rest-frame EW;
             # Note about 16/84 in fbb
-            EW16 = (fy_ex * c / np.square(x_ex) / d - fbb84[iix2]) / (fbb[iix2]) * lfwhm[iix2] / (1.+zbes)
-            EW50 = (fy_ex * c / np.square(x_ex) / d - fbb[iix2]) / (fbb[iix2]) * lfwhm[iix2] / (1.+zbes)
-            EW84 = (fy_ex * c / np.square(x_ex) / d - fbb16[iix2]) / (fbb[iix2]) * lfwhm[iix2] / (1.+zbes)
+            EW16 = (fy_ex * c / np.square(x_ex) /d_scale - fbb84[iix2]) / (fbb[iix2]) * lfwhm[iix2] / (1.+zbes)
+            EW50 = (fy_ex * c / np.square(x_ex) /d_scale - fbb[iix2]) / (fbb[iix2]) * lfwhm[iix2] / (1.+zbes)
+            EW84 = (fy_ex * c / np.square(x_ex) /d_scale - fbb16[iix2]) / (fbb[iix2]) * lfwhm[iix2] / (1.+zbes)
 
-            EW50_er1 = ((fy_ex-ey_ex) * c / np.square(x_ex) / d - fbb[iix2]) / (fbb[iix2]) * lfwhm[iix2] / (1.+zbes)
-            EW50_er2 = ((fy_ex+ey_ex) * c / np.square(x_ex) / d - fbb[iix2]) / (fbb[iix2]) * lfwhm[iix2] / (1.+zbes)
+            EW50_er1 = ((fy_ex-ey_ex) * c / np.square(x_ex) /d_scale - fbb[iix2]) / (fbb[iix2]) * lfwhm[iix2] / (1.+zbes)
+            EW50_er2 = ((fy_ex+ey_ex) * c / np.square(x_ex) /d_scale - fbb[iix2]) / (fbb[iix2]) * lfwhm[iix2] / (1.+zbes)
 
             cnt50 = fbb[iix2]
             cnt16 = fbb16[iix2]
@@ -2349,9 +2332,9 @@ def func_tmp(xint,eobs,fmodel):
 
         col5  = fits.Column(name='wave_obs', format='E', unit='AA', array=xbb)
         col00.append(col5)
-        col6  = fits.Column(name='f_obs', format='E', unit='1e%derg/s/cm2/AA'%(np.log10(scale)), array=fybb[:] * c / np.square(xbb[:]) / d)
+        col6  = fits.Column(name='f_obs', format='E', unit='1e%derg/s/cm2/AA'%(np.log10(scale)), array=fybb[:] * c / np.square(xbb[:]) /d_scale)
         col00.append(col6)
-        col7  = fits.Column(name='e_obs', format='E', unit='1e%derg/s/cm2/AA'%(np.log10(scale)), array=eybb[:] * c / np.square(xbb[:]) / d)
+        col7  = fits.Column(name='e_obs', format='E', unit='1e%derg/s/cm2/AA'%(np.log10(scale)), array=eybb[:] * c / np.square(xbb[:]) /d_scale)
         col00.append(col7)
 
         hdr = fits.Header()
@@ -2500,18 +2483,18 @@ def func_tmp(xint,eobs,fmodel):
             tree_spec.update({'f_model_dust': ytmp_dust50})            
         # BB for dust
         tree_spec.update({'wave_obs': xbb})
-        tree_spec.update({'f_obs': fybb[:] * c / np.square(xbb[:]) / d})
-        tree_spec.update({'e_obs': eybb[:] * c / np.square(xbb[:]) / d})
+        tree_spec.update({'f_obs': fybb[:] * c / np.square(xbb[:]) /d_scale})
+        tree_spec.update({'e_obs': eybb[:] * c / np.square(xbb[:]) /d_scale})
         # grism:
         if f_grsm:
-            tree_spec.update({'fg0_obs': fg0 * c/np.square(xg0)/d})
-            tree_spec.update({'eg0_obs': eg0 * c/np.square(xg0)/d})
+            tree_spec.update({'fg0_obs': fg0 * c/np.square(xg0)/d_scale})
+            tree_spec.update({'eg0_obs': eg0 * c/np.square(xg0)/d_scale})
             tree_spec.update({'wg0_obs': xg0})
-            tree_spec.update({'fg1_obs': fg1 * c/np.square(xg1)/d})
-            tree_spec.update({'eg1_obs': eg1 * c/np.square(xg1)/d})
+            tree_spec.update({'fg1_obs': fg1 * c/np.square(xg1)/d_scale})
+            tree_spec.update({'eg1_obs': eg1 * c/np.square(xg1)/d_scale})
             tree_spec.update({'wg1_obs': xg1})
-            tree_spec.update({'fg2_obs': fg2 * c/np.square(xg2)/d})
-            tree_spec.update({'eg2_obs': eg2 * c/np.square(xg2)/d})
+            tree_spec.update({'fg2_obs': fg2 * c/np.square(xg2)/d_scale})
+            tree_spec.update({'eg2_obs': eg2 * c/np.square(xg2)/d_scale})
             tree_spec.update({'wg2_obs': xg2})
 
         af = asdf.AsdfFile(tree_spec)
@@ -2543,8 +2526,8 @@ def func_tmp(xint,eobs,fmodel):
         ax2t.set_xlim(xgmin, xgmax)
 
         conaa = (x0>xgmin-50) & (x0<xgmax+50)
-        ymaxzoom = np.max(ysum[conaa]*c/np.square(x0[conaa])/d) * 1.15
-        yminzoom = np.min(ysum[conaa]*c/np.square(x0[conaa])/d) / 1.15
+        ymaxzoom = np.max(ysum[conaa]*c/np.square(x0[conaa])/d_scale) * 1.15
+        yminzoom = np.min(ysum[conaa]*c/np.square(x0[conaa])/d_scale) / 1.15
 
         ax2t.set_ylim(yminzoom, ymaxzoom)
         ax2t.xaxis.labelpad = -2
@@ -2558,7 +2541,7 @@ def func_tmp(xint,eobs,fmodel):
     if f_dust:
         try:
             contmp = (x1_tot>10*1e4)
-            y3min, y3max = -.2*np.max((model_tot * c/ np.square(x1_tot) / d)[contmp]), np.max((model_tot * c/ np.square(x1_tot) / d)[contmp])*2.0
+            y3min, y3max = -.2*np.max((model_tot * c/ np.square(x1_tot) /d_scale)[contmp]), np.max((model_tot * c/ np.square(x1_tot) /d_scale)[contmp])*2.0
             ax3t.set_ylim(y3min, y3max)
         except:
             if verbose:
@@ -2672,7 +2655,7 @@ def plot_filter(MB, ax, ymax, scl=0.3, cmap='gist_rainbow', alp=0.4, ind_remove=
 
 
 def plot_corner_physparam_summary(MB, fig=None, out_ind=0, DIR_OUT='./', mmax:int=1000, TMIN=0.0001, tau_lim=0.01, f_plot_filter=True, 
-    scale=1e-19, NRbb_lim=10000, save_pcl=True, return_figure=False):
+    scale=1e-19, NRbb_lim=10000, save_pcl=True, return_figure=False, SNlim=1):
     '''
     Purpose
     -------
@@ -2836,20 +2819,21 @@ def plot_corner_physparam_summary(MB, fig=None, out_ind=0, DIR_OUT='./', mmax:in
     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_scale = MB.d * scale
+    d = d_scale
 
     # BB photometry
     conspec = (NR<NRbb_lim)
     sigma = 1.
     conbb = (fybb/eybb > sigma)
-    ax0.errorbar(xbb[conbb], fybb[conbb] * c / np.square(xbb[conbb]) / d, yerr=eybb[conbb]*c/np.square(xbb[conbb])/d, color='k', linestyle='', linewidth=0.5, zorder=4)
-    ax0.plot(xbb[conbb], fybb[conbb] * c / np.square(xbb[conbb]) / d, '.r', ms=10, linestyle='', linewidth=0, zorder=4)
+    ax0.errorbar(xbb[conbb], fybb[conbb] * c / np.square(xbb[conbb]) /d_scale, yerr=eybb[conbb]*c/np.square(xbb[conbb])/d_scale, color='k', linestyle='', linewidth=0.5, zorder=4)
+    ax0.plot(xbb[conbb], fybb[conbb] * c / np.square(xbb[conbb]) /d_scale, '.r', ms=10, linestyle='', linewidth=0, zorder=4)
 
     conebb_ls = (fybb/eybb <= sigma)
     if len(fybb[conebb_ls])>0:
-        leng = np.max(fybb[conebb_ls] * c / np.square(xbb[conebb_ls]) / d) * 0.05
-        ax0.errorbar(xbb[conebb_ls], eybb[conebb_ls] * c / np.square(xbb[conebb_ls]) / d * sigma, yerr=leng,\
-            uplims=eybb[conebb_ls] * c / np.square(xbb[conebb_ls]) / d * sigma, linestyle='',color='r', marker='', ms=4, label='', zorder=4, capsize=3)
+        leng = np.max(fybb[conebb_ls] * c / np.square(xbb[conebb_ls]) /d_scale) * 0.05
+        ax0.errorbar(xbb[conebb_ls], eybb[conebb_ls] * c / np.square(xbb[conebb_ls]) /d_scale * sigma, yerr=leng,\
+            uplims=eybb[conebb_ls] * c / np.square(xbb[conebb_ls]) /d_scale * sigma, linestyle='',color='r', marker='', ms=4, label='', zorder=4, capsize=3)
     
     ####################
     # MCMC corner plot.
@@ -3104,14 +3088,14 @@ def density_estimation(m1, m2):
                 ysump = y0p
                 ysum = y0
                 if AA_tmp/Asum > flim:
-                    ax0.plot(x0, y0 * c/ np.square(x0) / d, '--', lw=.1, color=col[ii], zorder=-1, label='', alpha=0.01)
+                    ax0.plot(x0, y0 * c/ np.square(x0) /d_scale, '--', lw=.1, color=col[ii], zorder=-1, label='', alpha=0.01)
             else:
                 y0_r, x0_tmp = fnc.get_template_single(AA_tmp, Avtmp[kk], ii, ZZ_tmp, zbes, lib_all)
                 y0p, x0p = fnc.get_template_single(AA_tmp, Avtmp[kk], ii, ZZ_tmp, zbes, lib)
                 ysump += y0p
                 ysum += y0_r
                 if AA_tmp/Asum > flim:
-                    ax0.plot(x0, y0_r * c/ np.square(x0) / d, '--', lw=.1, color=col[ii], zorder=-1, label='', alpha=0.01)
+                    ax0.plot(x0, y0_r * c/ np.square(x0) /d_scale, '--', lw=.1, color=col[ii], zorder=-1, label='', alpha=0.01)
 
         for ss in range(len(age)):
             ii = ss # from old to young templates.
@@ -3122,7 +3106,7 @@ def density_estimation(m1, m2):
                 ZC[ss] = -99
 
         # Plot Total
-        ax0.plot(x0, ysum * c/ np.square(x0) / d, '-', lw=.1, color='gray', zorder=-1, label='', alpha=0.01)
+        ax0.plot(x0, ysum * c/ np.square(x0) /d_scale, '-', lw=.1, color='gray', zorder=-1, label='', alpha=0.01)
 
         if len(age)==1:
             ax1.plot(age[:], SF[:], marker='.', linestyle='-', lw=.1, color='k', zorder=-1, label='', alpha=0.01)
@@ -3137,8 +3121,8 @@ def density_estimation(m1, m2):
         else:
             scl_yaxis = 0
 
-        ymax_bb = np.max(fybb[conbb] * c / np.square(xbb[conbb]) / d) * 1.10
-        ymax_temp = np.max(ysum * c/ np.square(x0) / d) * 1.10
+        ymax_bb = np.max(fybb[conbb] * c / np.square(xbb[conbb]) /d_scale) * 1.10
+        ymax_temp = np.max(ysum * c/ np.square(x0) /d_scale) * 1.10
         ymax = np.max([ymax_bb, ymax_temp])
 
         # Convert into log
@@ -3260,9 +3244,9 @@ def density_estimation(m1, m2):
             files.append(fname)
 
     # For the last one
-    ax0.plot(xg0, fg0 * c / np.square(xg0) / d, marker='', linestyle='-', linewidth=0.5, ms=0.1, color='royalblue', label='')
-    ax0.plot(xg1, fg1 * c / np.square(xg1) / d, marker='', linestyle='-', linewidth=0.5, ms=0.1, color='g', label='')
-    ax0.plot(xg2, fg2 * c / np.square(xg2) / d, marker='', linestyle='-', linewidth=0.5, ms=0.1, color='#DF4E00', label='')
+    ax0.plot(xg0, fg0 * c / np.square(xg0) /d_scale, marker='', linestyle='-', linewidth=0.5, ms=0.1, color='royalblue', label='')
+    ax0.plot(xg1, fg1 * c / np.square(xg1) /d_scale, marker='', linestyle='-', linewidth=0.5, ms=0.1, color='g', label='')
+    ax0.plot(xg2, fg2 * c / np.square(xg2) /d_scale, marker='', linestyle='-', linewidth=0.5, ms=0.1, color='#DF4E00', label='')
 
     ax0.set_xlim(2200, 88000)
     ax0.set_xscale('log')
@@ -3369,8 +3353,8 @@ def write_lines(ID, zbes, R_grs=45, dw=4, umag=1.0, ldw = 7, DIR_OUT='./'):
                 flux = np.zeros(len(xxs), dtype=float)
                 efl  = np.zeros(len(xxs), dtype=float)
                 for ff in range(len(xxs)):
-                    flux[ff] = yy2s[ff]/np.square(xxs[ff]) * c/d
-                    efl[ff]  = np.square(eyys[ff]/np.square(xxs[ff]) * c/d)
+                    flux[ff] = yy2s[ff]/np.square(xxs[ff]) * c/d_scale
+                    efl[ff]  = np.square(eyys[ff]/np.square(xxs[ff]) * c/d_scale)
 
                 fmed = np.median(flux) # Median of continuum, model flux
                 esum = np.sqrt(simps(efl, xxs))
@@ -3384,33 +3368,33 @@ def write_lines(ID, zbes, R_grs=45, dw=4, umag=1.0, ldw = 7, DIR_OUT='./'):
                         xxss = xxs/zscl
 
                     if f_grsm:
-                        ax2t.plot(xxs/zscl, (gaus(xxs,*popt)+yy2s) * c/np.square(xxs)/d, '#4682b4', linestyle='-', linewidth=1, alpha=0.8, zorder=20)
+                        ax2t.plot(xxs/zscl, (gaus(xxs,*popt)+yy2s) * c/np.square(xxs)/d_scale, '#4682b4', linestyle='-', linewidth=1, alpha=0.8, zorder=20)
 
-                    I1 = simps((gaus(xxs,*popt)) * c/np.square(xxs)/d, xxs)
-                    I2 = I1 - simps((gaus(xxs,*popt)) * c/np.square(xxs)/d, xxs)
+                    I1 = simps((gaus(xxs,*popt)) * c/np.square(xxs)/d_scale, xxs)
+                    I2 = I1 - simps((gaus(xxs,*popt)) * c/np.square(xxs)/d_scale, xxs)
                     fline = I1
 
                     Flum = fline*Cons*1e-18 # luminosity in erg/s.
                     elum = esum *Cons*1e-18 # luminosity in erg/s.
                     SFR  = Flum * 6.58*1e-42
                     print('SFR is', SFR/umag)
-                    EW_tmp   = simps( ((gaus(xxs,*popt)) * c/np.square(xxs)/d)/yy2s, xxs)
-                    EW_tmp_u = simps( ((gaus(xxs,*popt) + eyys/np.sqrt(len(xxs))) * c/np.square(xxs)/d)/yy2s, xxs)
+                    EW_tmp   = simps( ((gaus(xxs,*popt)) * c/np.square(xxs)/d_scale)/yy2s, xxs)
+                    EW_tmp_u = simps( ((gaus(xxs,*popt) + eyys/np.sqrt(len(xxs))) * c/np.square(xxs)/d_scale)/yy2s, xxs)
 
                     if ii == 7:
                         contmp2 = (xxs/zscl>4320.) & (xxs/zscl<4380.)
                         popt,pcov = curve_fit(gaus,xxs[contmp2], yys[contmp2], p0=[Fline50[ii],WL,10], sigma=eyys[contmp2])
 
-                        I1 = simps((gaus(xxs[contmp2],*popt)) * c/np.square(xxs[contmp2])/d, xxs[contmp2])
-                        I2 = I1 - simps((gaus(xxs[contmp2],*popt)) * c/np.square(xxs[contmp2])/d, xxs[contmp2])
+                        I1 = simps((gaus(xxs[contmp2],*popt)) * c/np.square(xxs[contmp2])/d_scale, xxs[contmp2])
+                        I2 = I1 - simps((gaus(xxs[contmp2],*popt)) * c/np.square(xxs[contmp2])/d_scale, xxs[contmp2])
                         fline = I1
 
                         Flum = fline*Cons*1e-18 # luminosity in erg/s.
                         elum = esum *Cons*1e-18 # luminosity in erg/s.
                         SFR  = Flum * 6.58*1e-42
                         print('SFR, update, is', SFR/umag)
-                        EW_tmp   = simps( ((gaus(xxs[contmp2],*popt)) * c/np.square(xxs[contmp2])/d)/yy2s[contmp2], xxs[contmp2])
-                        EW_tmp_u = simps( ((gaus(xxs[contmp2],*popt) + eyys[contmp2]/np.sqrt(len(xxs[contmp2]))) * c/np.square(xxs[contmp2])/d)/yy2s[contmp2], xxs[contmp2])
+                        EW_tmp   = simps( ((gaus(xxs[contmp2],*popt)) * c/np.square(xxs[contmp2])/d_scale)/yy2s[contmp2], xxs[contmp2])
+                        EW_tmp_u = simps( ((gaus(xxs[contmp2],*popt) + eyys[contmp2]/np.sqrt(len(xxs[contmp2]))) * c/np.square(xxs[contmp2])/d_scale)/yy2s[contmp2], xxs[contmp2])
 
                     flw.write('%d %.2f %.2f %.2f %.2f %.2f %.2e %.2e %.2f\n'%(LW[ii],fline/umag, esum/umag, fmed/umag, EW_tmp,(EW_tmp_u-EW_tmp), Flum*1e-18/umag, elum*1e-18/umag, SFR/umag))
 
@@ -3419,7 +3403,7 @@ def write_lines(ID, zbes, R_grs=45, dw=4, umag=1.0, ldw = 7, DIR_OUT='./'):
                     for ff in range(len(fsum)):
                         fsum[ff] = (yys[ff]+yy2s[ff])/np.square(xxs[ff])
 
-                    fline = np.sum(fsum) / d*c
+                    fline = np.sum(fsum) /d_scale*c
                     flw.write('%d %.2f %.2f %.2f %d %d %d %d %d\n'%(LW[ii],fline,esum,fmed, -99, 0, -99, 0, 0))
                     pass
 
@@ -3551,8 +3535,8 @@ def plot_corner_physparam_frame(ID, PA, Zall=np.arange(-1.2,0.4249,0.05), age=[0
     snbb = fybb/eybb
 
     conspec = (NR<10000) #& (fy/ey>1)
-    #ax0.plot(xg0, fg0 * c / np.square(xg0) / d, marker='', linestyle='-', linewidth=0.5, ms=0.1, color='royalblue', label='')
-    #ax0.plot(xg1, fg1 * c / np.square(xg1) / d, marker='', linestyle='-', linewidth=0.5, ms=0.1, color='#DF4E00', label='')
+    #ax0.plot(xg0, fg0 * c / np.square(xg0) /d_scale, marker='', linestyle='-', linewidth=0.5, ms=0.1, color='royalblue', label='')
+    #ax0.plot(xg1, fg1 * c / np.square(xg1) /d_scale, marker='', linestyle='-', linewidth=0.5, ms=0.1, color='#DF4E00', label='')
     conbb = (fybb/eybb>1)
 
     DIR_TMP = './templates/'
@@ -3675,10 +3659,10 @@ def plot_corner_physparam_frame(ID, PA, Zall=np.arange(-1.2,0.4249,0.05), age=[0
             ax1 = fig.add_axes([0.05,0.40,0.37,0.23])
             ax2 = fig.add_axes([0.05,0.07,0.37,0.23])
 
-            ax0.errorbar(xbb[conbb], fybb[conbb] * c / np.square(xbb[conbb]) / d, yerr=eybb[conbb]*c/np.square(xbb[conbb])/d, color='k', linestyle='', linewidth=0.5, zorder=4)
-            ax0.plot(xbb[conbb], fybb[conbb] * c / np.square(xbb[conbb]) / d, '.r', ms=10, linestyle='', linewidth=0, zorder=4)#, label='Obs.(BB)')
-            ax0.plot(xg0, fg0 * c / np.square(xg0) / d, marker='', linestyle='-', linewidth=0.5, ms=0.1, color='royalblue', label='')
-            ax0.plot(xg1, fg1 * c / np.square(xg1) / d, marker='', linestyle='-', linewidth=0.5, ms=0.1, color='#DF4E00', label='')
+            ax0.errorbar(xbb[conbb], fybb[conbb] * c / np.square(xbb[conbb]) /d_scale, yerr=eybb[conbb]*c/np.square(xbb[conbb])/d_scale, color='k', linestyle='', linewidth=0.5, zorder=4)
+            ax0.plot(xbb[conbb], fybb[conbb] * c / np.square(xbb[conbb]) /d_scale, '.r', ms=10, linestyle='', linewidth=0, zorder=4)#, label='Obs.(BB)')
+            ax0.plot(xg0, fg0 * c / np.square(xg0) /d_scale, marker='', linestyle='-', linewidth=0.5, ms=0.1, color='royalblue', label='')
+            ax0.plot(xg1, fg1 * c / np.square(xg1) /d_scale, marker='', linestyle='-', linewidth=0.5, ms=0.1, color='#DF4E00', label='')
 
 
         #nr = kk # np.random.randint(len(samples))
@@ -3723,7 +3707,7 @@ def plot_corner_physparam_frame(ID, PA, Zall=np.arange(-1.2,0.4249,0.05), age=[0
                 ysump = y0p #* 1e18
                 ysum  = y0  #* 1e18
                 if AA_tmp/Asum > flim:
-                    ax0.plot(x0, y0 * c/ np.square(x0) / d, '--', lw=1, color=col[ii], zorder=-1, label='', alpha=0.5)
+                    ax0.plot(x0, y0 * c/ np.square(x0) /d_scale, '--', lw=1, color=col[ii], zorder=-1, label='', alpha=0.5)
                 #ax1.plot(age[ii], SF[ii], marker='o', lw=1, color=col[ii], zorder=1, label='', alpha=0.5)
                 #ax1.errorbar(age[ii], SF[ii], xerr=[[delTl[ii]/1e9], [delTu[ii]/1e9]], ms=10, marker='o', lw=1, color=col[ii], zorder=1, label='', alpha=0.5)
                 xx1 = np.arange(age[ii]-delTl[ii]/1e9, age[ii]+delTu[ii]/1e9, 0.01)
@@ -3735,7 +3719,7 @@ def plot_corner_physparam_frame(ID, PA, Zall=np.arange(-1.2,0.4249,0.05), age=[0
                 ysump += y0p  #* 1e18
                 ysum  += y0_r #* 1e18
                 if AA_tmp/Asum > flim:
-                    ax0.plot(x0, y0_r * c/ np.square(x0) / d, '--', lw=1, color=col[ii], zorder=-1, label='', alpha=0.5)
+                    ax0.plot(x0, y0_r * c/ np.square(x0) /d_scale, '--', lw=1, color=col[ii], zorder=-1, label='', alpha=0.5)
                 #ax1.plot(age[ii], SF[ii], marker='o', lw=1, color=col[ii], zorder=1, label='', alpha=0.5)
                 #ax1.errorbar(age[ii], SF[ii], xerr=[[delTl[ii]/1e9], [delTu[ii]/1e9]], ms=10, marker='o', lw=1, color=col[ii], zorder=1, label='', alpha=0.5)
                 xx1 = np.arange(age[ii]-delTl[ii]/1e9, age[ii]+delTu[ii]/1e9, 0.01)
@@ -3751,8 +3735,8 @@ def plot_corner_physparam_frame(ID, PA, Zall=np.arange(-1.2,0.4249,0.05), age=[0
             ax2.errorbar(age[ii], ZC[ii], xerr=[[delTl[ii]/1e9], [delTu[ii]/1e9]], ms=10*(AC/np.sum(AM[:]))+1, marker='o', lw=1, color=col[ii], zorder=1, label='', alpha=0.5)
 
         # Total
-        ymax = np.max(fybb[conbb] * c / np.square(xbb[conbb]) / d) * 1.10
-        ax0.plot(x0, ysum * c/ np.square(x0) / d, '-', lw=1., color='gray', zorder=-1, label='', alpha=0.8)
+        ymax = np.max(fybb[conbb] * c / np.square(xbb[conbb]) /d_scale) * 1.10
+        ax0.plot(x0, ysum * c/ np.square(x0) /d_scale, '-', lw=1., color='gray', zorder=-1, label='', alpha=0.8)
         ax0.set_xlim(2200, 88000)
         ax0.set_xscale('log')
         ax0.set_ylim(0., ymax)
diff --git a/gsf/plot_sfh.py b/gsf/plot_sfh.py
index 6d6f1a7..e3f4f8e 100644
--- a/gsf/plot_sfh.py
+++ b/gsf/plot_sfh.py
@@ -42,7 +42,8 @@ def plot_sfh(MB, flim=0.01, lsfrl=-3, mmax=1000, Txmin=0.08, Txmax=4, lmmin=5, f
         import matplotlib
         matplotlib.use("Agg")
 
-    print('\n### Running plot_sfh ###\n')
+    MB.logger.info('Running plot_sfh')
+
     fnc = MB.fnc
     bfnc = MB.bfnc
     ID = MB.ID
@@ -276,7 +277,7 @@ def plot_sfh(MB, flim=0.01, lsfrl=-3, mmax=1000, Txmin=0.08, Txmax=4, lmmin=5, f
                 pass
     except:
         if verbose:
-            print('No simulation file (%s).\nError may be underestimated.' % meanfile)
+            MB.logger.warning('No simulation file (%s).\nError may be underestimated.' % meanfile)
         eA = age * 0
         eZ = age * 0
         eAv = 0
@@ -470,7 +471,7 @@ def plot_sfh(MB, flim=0.01, lsfrl=-3, mmax=1000, Txmin=0.08, Txmax=4, lmmin=5, f
         f_rejuv,t_quench,t_rejuv = check_rejuv(age,SFp[:,:],ACp[:,:],SFMS_50)
     else:
         if verbose:
-            print('Failed to call rejuvenation module.')
+            MB.logger.warning('Failed to call rejuvenation module.')
         f_rejuv,t_quench,t_rejuv = 0,0,0
 
     # Plot MS?
diff --git a/gsf/writing.py b/gsf/writing.py
index 0c93fc4..e9da20b 100644
--- a/gsf/writing.py
+++ b/gsf/writing.py
@@ -24,12 +24,12 @@ def get_param(self, res, fitc, tcalc=1., burnin=-1):
     Czrec1 = self.Cz1
     Czrec2 = self.Cz2
 
-    try:
+    if self.fzvis:
         z_cz = self.z_cz_prev
         scl_cz0 = self.scl_cz0
         scl_cz1 = self.scl_cz1
         scl_cz2 = self.scl_cz2
-    except: # When redshiftfit is skipped.
+    else: # When redshiftfit is skipped.
         z_cz = np.asarray([self.zgal,self.zgal,self.zgal])
         scl_cz0 = np.asarray([self.Cz0,self.Cz0,self.Cz0])
         scl_cz1 = np.asarray([self.Cz1,self.Cz1,self.Cz1])