Merge pull request #36 from mtakahiro/nirspec-fit

Full Spectral Fitting
mtakahiro · Nov 23, 2022 · ec9fb19 · ec9fb19
2 parents 61cae9a + 9cf5fef
commit ec9fb19
Show file tree

Hide file tree

Showing 11 changed files with 383 additions and 190 deletions.
diff --git a/docs/parameters.rst b/docs/parameters.rst
@@ -136,6 +136,9 @@ An example config file can be generated by executing the following command (TBD)
    * - AMIN
      - str 
      - (Optional) Minimum value for amplitude, in normal logarithmic scale.
+   * - DUST_MODEL
+     - int
+     - Dust attenuation model (default 0). 0 for Calzetti. 1 for MW. 2 for LCM. 3 for SMC. 4 for Kriek & Conroy (2013).
    * - AVMAX
      - float 
      - (Optional) Maximum value for Av (dust attenuation in V-band), in mag.
@@ -171,7 +174,7 @@ An example config file can be generated by executing the following command (TBD)
      - 
 
 
-**Parameters for a specific target:**
+**Parameters specific for the target:**
 
 .. list-table::
    :widths: 10 5 20

diff --git a/gsf/fitting.py b/gsf/fitting.py
diff --git a/gsf/function.py b/gsf/function.py
@@ -248,7 +248,8 @@ def loadcpkl(cpklfile):
     return data
 
 
-def get_leastsq(MB, ZZtmp, fneld, age, fit_params, residual, fy, ey, wht, ID0, chidef=None, Zbest=0, f_keep=False):
+def get_leastsq(MB, ZZtmp, fneld, age, fit_params, residual, fy, ey, wht, ID0, 
+    chidef=None, Zbest=0, f_keep=False, f_only_spec=False):
     '''
     Get initial parameters at various Z
     '''
@@ -283,7 +284,10 @@ def get_leastsq(MB, ZZtmp, fneld, age, fit_params, residual, fy, ey, wht, ID0, c
                 fit_params['Z'+str(aa)].value = ZZ
 
         f_fir = False
-        out_tmp = minimize(residual, fit_params, args=(fy, ey, wht, f_fir), method=fit_name) # nelder is the most efficient.
+
+        out_tmp = minimize(residual, fit_params, args=(fy, ey, wht, f_fir), 
+            method=fit_name, kws={'f_only_spec':f_only_spec})
+
         csq = out_tmp.chisqr
         rcsq = out_tmp.redchi
         fitc = [csq, rcsq] # Chi2, Reduced-chi2
@@ -550,12 +554,18 @@ def flamtonu(lam, flam, m0set=25.0):
     return fnu
 
 
-def fnutolam(lam, fnu, m0set=25.0):
+def fnutolam(lam, fnu, m0set=25.0, m0=-48.6):
     '''
     Converts from Fnu to Flam, from mag zeropoint of m0set (to -48.6).
-    
+
+    Parameters
+    ----------
+    m0set : float
+        current magzp.
+    m0 : float
+        target magzp. The default, -48.6, is for flam (erg/s/cm2/lambda).
     '''
-    Ctmp = lam**2/c * 10**((48.6+m0set)/2.5)
+    Ctmp = lam**2/c * 10**((m0set-m0)/2.5)
     flam = fnu / Ctmp
     return flam
 

diff --git a/gsf/function_class.py b/gsf/function_class.py
@@ -224,7 +224,7 @@ def open_spec_fits_dir(self, nage:int, nz:int, kk, Av00:float, zgal:float, A00:f
         elif self.dust_model == 4: # Kriek&Conroy with gamma=-0.2
             yyd, xxd, nrd = dust_kc(xx, yy, Av00, nr, Rv=4.05, gamma=-0.2)
         else:
-            print('No entry. Dust model is set to Calzetti')
+            print('No valid entry. Dust model is set to Calzetti')
             yyd, xxd, nrd = dust_calz(xx, yy, Av00, nr)
 
         xxd *= (1.+zgal)
@@ -940,7 +940,7 @@ def open_spec_fits_dir(self, nage, nz, kk, Av00, zgal, A00):
         elif self.dust_model == 4: # Kriek&Conroy with gamma=-0.2
             yyd, xxd, nrd = dust_kc(xx, yy, Av00, nr, Rv=4.05, gamma=-0.2)
         else:
-            print('No entry. Dust model is set to Calzetti')
+            print('No valid entry. Dust model is set to Calzetti')
             yyd, xxd, nrd = dust_calz(xx, yy, Av00, nr)
 
         xxd *= (1.+zgal)

diff --git a/gsf/maketmp_filt.py b/gsf/maketmp_filt.py
@@ -435,16 +435,23 @@ def maketemp(MB, ebblim=1e10, lamliml=0., lamlimu=50000., ncolbb=10000,
     # READ BB photometry from CAT_BB:
     #############################
     if CAT_BB:
+        key_id = 'id'
         fd0 = ascii.read(CAT_BB)
-        id0 = fd0['id'].astype('str')
+        try:
+            id0 = fd0[key_id].astype('str')
+        except:
+            key_id = 'ID'
+            id0 = fd0[key_id].astype('str')
+
         ii0 = np.where(id0[:]==ID)
         try:
             if len(ii0[0]) == 0:
                 print('Could not find the column for [ID: %s] in the input BB catalog! Exiting.'%(ID))
                 return False
-            id = fd0['id'][ii0]
+            id = fd0[key_id][ii0]
         except:
             print('Could not find the column for [ID: %s] in the input BB catalog! Exiting.'%(ID))
+            print(fd0)
             return False
 
         fbb = np.zeros(len(SFILT), dtype='float')
@@ -494,18 +501,25 @@ def maketemp(MB, ebblim=1e10, lamliml=0., lamlimu=50000., ncolbb=10000,
 
     # Dust flux;
     if MB.f_dust:
+        key_id = 'id'
         fdd = ascii.read(CAT_BB_DUST)
         try:
-            id0 = fdd['id'].astype('str')
+            try:
+                id0 = fdd[key_id].astype('str')
+            except:
+                key_id = 'ID'
+                id0 = fdd[key_id].astype('str')
+
             ii0 = np.where(id0[:]==ID)
             try:
-                id = fd0['id'][ii0]
+                id = fd0[key_id][ii0]
             except:
                 print('Could not find the column for [ID: %s] in the input BB catalog! Exiting.'%(ID))
                 return False
         except:
             return False
-        id = fdd['id']
+
+        id = fdd[key_id]
 
         fbb_d = np.zeros(len(DFILT), dtype='float')
         ebb_d = np.zeros(len(DFILT), dtype='float')
@@ -617,7 +631,7 @@ def maketemp(MB, ebblim=1e10, lamliml=0., lamlimu=50000., ncolbb=10000,
                     try:
                         spec_mul_nu_conv[ss,:] = convolve(spec_mul_nu[ss], LSF, boundary='extend')
                     except:
-                        print('Error. No convolution is happening...')
+                        #print('Error. No convolution is happening...')
                         spec_mul_nu_conv[ss,:] = spec_mul_nu[ss]
                         if zz==0 and ss==0:
                             print('Kernel is too small. No convolution.')

diff --git a/gsf/maketmp_z0.py b/gsf/maketmp_z0.py
@@ -9,7 +9,7 @@
 INDICES = ['G4300', 'Mgb', 'Fe5270', 'Fe5335', 'NaD', 'Hb', 'Fe4668', 'Fe5015', 'Fe5709', 'Fe5782', 'Mg1', 'Mg2', 'TiO1', 'TiO2']
 
 
-def make_tmp_z0(MB, lammin=100, lammax=160000, tau_lim=0.001, force_no_neb=False, Zforce=None):
+def make_tmp_z0(MB, lammin=100, lammax=160000, tau_lim=0.001, force_no_neb=False, Zforce=None, f_mp=True):
     '''
     This is for the preparation of default template, with FSPS, at z=0.
     Should be run before SED fitting.

diff --git a/gsf/plot_sed.py b/gsf/plot_sed.py
@@ -17,6 +17,7 @@
 
 col = ['violet', 'indigo', 'b', 'lightblue', 'lightgreen', 'g', 'orange', 'coral', 'r', 'darkred']#, 'k']
 
+
 def plot_sed(MB, flim=0.01, fil_path='./', scale=1e-19, f_chind=True, figpdf=False, save_sed=True, 
     mmax=300, dust_model=0, DIR_TMP='./templates/', f_label=False, f_bbbox=False, verbose=False, f_silence=True,
     f_fill=False, f_fancyplot=False, f_Alog=True, dpi=300, f_plot_filter=True, f_plot_resid=False, NRbb_lim=10000):
@@ -98,7 +99,6 @@ def gaus(x,a,x0,sigma):
     chimax = 1.
     m0set = MB.m0set
     Mpc_cm = MB.Mpc_cm
-    d = MB.d * scale 
 
     ##################
     # Fitting Results
@@ -289,14 +289,20 @@ def gaus(x,a,x0,sigma):
             AAAA
             BBBB
             """
-            fig,axes = plt.subplot_mosaic(mosaic=fig_mosaic, figsize=(5.5,4.2))
-            fig.subplots_adjust(top=0.98, bottom=0.16, left=0.1, right=0.99, hspace=0.15, wspace=0.25)
+            fig,axes = plt.subplot_mosaic(mosaic=fig_mosaic, figsize=(5.5,4.))
+            fig.subplots_adjust(top=0.98, bottom=0.16, left=0.08, right=0.99, hspace=0.15, wspace=0.25)
             ax1 = axes['A']
         else:
-            fig = plt.figure(figsize=(5.5,2.2))
-            fig.subplots_adjust(top=0.98, bottom=0.16, left=0.1, right=0.99, hspace=0.15, wspace=0.25)
+            fig = plt.figure(figsize=(5.5,2.))
+            fig.subplots_adjust(top=0.98, bottom=0.16, left=0.08, right=0.99, hspace=0.15, wspace=0.25)
             ax1 = fig.add_subplot(111)
 
+    # Determine scale here;
+    if scale == None:
+        conbb_hs = (fybb/eybb > SNlim)
+        scale = 10**(int(np.log10(np.nanmax(fybb[conbb_hs] * c / np.square(xbb[conbb_hs])) / MB.d))) / 10
+    d = MB.d * scale
+
     #######################################
     # D.Kelson like Box for BB photometry
     #######################################
@@ -495,8 +501,8 @@ def gaus(x,a,x0,sigma):
     xboxl = 17000
     xboxu = 28000
 
-    ax1.set_xlabel('Observed wavelength ($\mathrm{\mu m}$)', fontsize=12)
-    ax1.set_ylabel('Flux ($10^{%d}\mathrm{erg}/\mathrm{s}/\mathrm{cm}^{2}/\mathrm{\AA}$)'%(np.log10(scale)),fontsize=12,labelpad=-2)
+    ax1.set_xlabel('Observed wavelength ($\mathrm{\mu m}$)', fontsize=11)
+    ax1.set_ylabel('Flux ($10^{%d}\mathrm{erg}/\mathrm{s}/\mathrm{cm}^{2}/\mathrm{\AA}$)'%(np.log10(scale)),fontsize=11,labelpad=-2)
 
     x1min = 2000
     x1max = 100000
@@ -778,17 +784,15 @@ def gaus(x,a,x0,sigma):
         from astropy.convolution import convolve
         from .maketmp_filt import get_LSF
         LSF, _ = get_LSF(MB.inputs, MB.DIR_EXTR, ID, x1_tot[:]/(1.+zbes), c=3e18)
-        spec_grsm16 = convolve(ytmp16[:], LSF, boundary='extend')
-        spec_grsm50 = convolve(ytmp50[:], LSF, boundary='extend')
-        spec_grsm84 = convolve(ytmp84[:], LSF, boundary='extend')
-        '''
-        plt.close()
-        plt.plot(x1_tot, ytmp50, color='r')
-        plt.plot(x1_tot, spec_grsm50, color='gray')
-        print(spec_grsm50)
-        plt.xlim(8000,20000)
-        plt.show()
-        '''
+        try:
+            spec_grsm16 = convolve(ytmp16[:], LSF, boundary='extend')
+            spec_grsm50 = convolve(ytmp50[:], LSF, boundary='extend')
+            spec_grsm84 = convolve(ytmp84[:], LSF, boundary='extend')
+        except:
+            spec_grsm16 = ytmp16[:]
+            spec_grsm50 = ytmp50[:]
+            spec_grsm84 = ytmp84[:]
+
         if True:
             ax2t.plot(x1_tot[:], ytmp50, '-', lw=0.5, color='gray', zorder=3., alpha=1.0)
         else:
@@ -1012,11 +1016,11 @@ def func_tmp(xint,eobs,fmodel):
 
         # Grism;
         if f_grsm:
-            col2  = fits.Column(name='f_model_conv_16', format='E', unit='1e%derg/s/cm2/AA'%(np.log10(scale)), array=spec_grsm16)
+            col2 = fits.Column(name='f_model_conv_16', format='E', unit='1e%derg/s/cm2/AA'%(np.log10(scale)), array=spec_grsm16)
             col00.append(col2)
-            col3  = fits.Column(name='f_model_conv_50', format='E', unit='1e%derg/s/cm2/AA'%(np.log10(scale)), array=spec_grsm50)
+            col3 = fits.Column(name='f_model_conv_50', format='E', unit='1e%derg/s/cm2/AA'%(np.log10(scale)), array=spec_grsm50)
             col00.append(col3)
-            col4  = fits.Column(name='f_model_conv_84', format='E', unit='1e%derg/s/cm2/AA'%(np.log10(scale)), array=spec_grsm84)
+            col4 = fits.Column(name='f_model_conv_84', format='E', unit='1e%derg/s/cm2/AA'%(np.log10(scale)), array=spec_grsm84)
             col00.append(col4)
 
         # BB for dust
@@ -1055,10 +1059,10 @@ def func_tmp(xint,eobs,fmodel):
 
         try:
             # Muv
-            MUV = -2.5 * np.log10(Fuv[:]) + 25.0
-            hdr['MUV16'] = -2.5 * np.log10(np.percentile(Fuv[:],16)) + 25.0
-            hdr['MUV50'] = -2.5 * np.log10(np.percentile(Fuv[:],50)) + 25.0
-            hdr['MUV84'] = -2.5 * np.log10(np.percentile(Fuv[:],84)) + 25.0
+            MUV = -2.5 * np.log10(Fuv[:]) + MB.m0set
+            hdr['MUV16'] = -2.5 * np.log10(np.percentile(Fuv[:],16)) + MB.m0set
+            hdr['MUV50'] = -2.5 * np.log10(np.percentile(Fuv[:],50)) + MB.m0set
+            hdr['MUV84'] = -2.5 * np.log10(np.percentile(Fuv[:],84)) + MB.m0set
 
             # Fuv (!= flux of Muv)
             hdr['FUV16'] = np.percentile(Fuv28[:],16)
@@ -1198,12 +1202,12 @@ def func_tmp(xint,eobs,fmodel):
     if f_label:
         fd = fits.open(MB.DIR_OUT + 'SFH_' + ID + '.fits')[0].header
         if MB.f_dust:
-            label = 'ID: %s\n$z:%.2f$\n$\log M_\mathrm{*}/M_\odot:%.2f$\n$\log M_\mathrm{dust}/M_\odot:%.2f$\n$T_\mathrm{dust}/K:%.1f$\n$\log Z_\mathrm{*}/Z_\odot:%.2f$\n$\log T_\mathrm{*}$/Gyr$:%.2f$\n$A_V$/mag$:%.2f$\n$\\chi^2/\\nu:%.2f$'\
-            %(ID, zbes, float(fd['Mstel_50']), MD50, TD50, float(fd['Z_MW_50']), float(fd['T_MW_50']), float(fd['AV_50']), fin_chi2)
+            label = 'ID: %s\n$z:%.2f$\n$\log M_\mathrm{*}/M_\odot:%.2f$\n$\log M_\mathrm{dust}/M_\odot:%.2f$\n$T_\mathrm{dust}/K:%.1f$\n$\log Z_\mathrm{*}/Z_\odot:%.2f$\n$\log T_\mathrm{*}$/Gyr$:%.2f$\n$A_V$/mag$:%.2f$'\
+            %(ID, zbes, float(fd['Mstel_50']), MD50, TD50, float(fd['Z_MW_50']), float(fd['T_MW_50']), float(fd['AV_50']))#, fin_chi2)
             ylabel = ymax*0.45
         else:
-            label = 'ID: %s\n$z:%.2f$\n$\log M_\mathrm{*}/M_\odot:%.2f$\n$\log Z_\mathrm{*}/Z_\odot:%.2f$\n$\log T_\mathrm{*}$/Gyr$:%.2f$\n$A_V$/mag$:%.2f$\n$\\chi^2/\\nu:%.2f$'\
-            %(ID, zbes, float(fd['Mstel_50']), float(fd['Z_MW_50']), float(fd['T_MW_50']), float(fd['AV_50']), fin_chi2)
+            label = 'ID: %s\n$z:%.2f$\n$\log M_\mathrm{*}/M_\odot:%.2f$\n$\log Z_\mathrm{*}/Z_\odot:%.2f$\n$\log T_\mathrm{*}$/Gyr$:%.2f$\n$A_V$/mag$:%.2f$'\
+            %(ID, zbes, float(fd['Mstel_50']), float(fd['Z_MW_50']), float(fd['T_MW_50']), float(fd['AV_50']))
             ylabel = ymax*0.25
 
         if f_grsm:
@@ -1402,7 +1406,6 @@ def gaus(x,a,x0,sigma):
     chimax = 1.
     m0set = MB.m0set
     Mpc_cm = MB.Mpc_cm
-    d = MB.d * scale
 
     ##################
     # Fitting Results
@@ -1614,6 +1617,12 @@ def gaus(x,a,x0,sigma):
             fig.subplots_adjust(top=0.98, bottom=0.16, left=0.1, right=0.99, hspace=0.15, wspace=0.25)
             ax1 = fig.add_subplot(111)
 
+    # Determine scale here;
+    if scale == None:
+        conbb_hs = (fybb/eybb > SNlim)
+        scale = 10**(int(np.log10(np.nanmax(fybb[conbb_hs] * c / np.square(xbb[conbb_hs])) / MB.d)))
+    d = MB.d * scale
+
     #######################################
     # D.Kelson like Box for BB photometry
     #######################################
@@ -2311,10 +2320,10 @@ def func_tmp(xint,eobs,fmodel):
 
         try:
             # Muv
-            MUV = -2.5 * np.log10(Fuv[:]) + 25.0
-            hdr['MUV16'] = -2.5 * np.log10(np.percentile(Fuv[:],16)) + 25.0
-            hdr['MUV50'] = -2.5 * np.log10(np.percentile(Fuv[:],50)) + 25.0
-            hdr['MUV84'] = -2.5 * np.log10(np.percentile(Fuv[:],84)) + 25.0
+            MUV = -2.5 * np.log10(Fuv[:]) + MB.m0set
+            hdr['MUV16'] = -2.5 * np.log10(np.percentile(Fuv[:],16)) + MB.m0set
+            hdr['MUV50'] = -2.5 * np.log10(np.percentile(Fuv[:],50)) + MB.m0set
+            hdr['MUV84'] = -2.5 * np.log10(np.percentile(Fuv[:],84)) + MB.m0set
 
             # Fuv (!= flux of Muv)
             hdr['FUV16'] = np.percentile(Fuv28[:],16)
@@ -2621,7 +2630,6 @@ def plot_corner_physparam_summary(MB, fig=None, out_ind=0, DIR_OUT='./', mmax=10
     ID   = MB.ID
     Z    = MB.Zall
     age  = MB.age
-    d = MB.d * scale
     c = MB.c
 
     tau0 = MB.tau0 
@@ -2762,6 +2770,12 @@ def plot_corner_physparam_summary(MB, fig=None, out_ind=0, DIR_OUT='./', mmax=10
     ey = np.append(ey01,eg_bb)
     wht = 1./np.square(ey)
 
+    # Determine scale here;
+    if scale == None:
+        conbb_hs = (fybb/eybb > SNlim)
+        scale = 10**(int(np.log10(np.nanmax(fybb[conbb_hs] * c / np.square(xbb[conbb_hs])) / MB.d)))
+    d = MB.d * scale
+
     # BB photometry
     conspec = (NR<NRbb_lim)
     sigma = 1.