diff --git a/doc/changes.rst b/doc/changes.rst index 38a44cec..dbab6315 100644 --- a/doc/changes.rst +++ b/doc/changes.rst @@ -23,6 +23,8 @@ Change Log * Switch from (deprecated) ``pkg_resources`` to ``importlib``. * Updated documentation (data model) and several non-negligible speed-ups. +* Improved modeling of galaxies with broad+narrow line-emission [`PR #142`_]: + .. _`PR #115`: https://github.com/desihub/fastspecfit/pull/115 .. _`PR #116`: https://github.com/desihub/fastspecfit/pull/116 .. _`PR #120`: https://github.com/desihub/fastspecfit/pull/120 @@ -30,6 +32,7 @@ Change Log .. _`PR #135`: https://github.com/desihub/fastspecfit/pull/135 .. _`PR #136`: https://github.com/desihub/fastspecfit/pull/136 .. _`PR #137`: https://github.com/desihub/fastspecfit/pull/137 +.. _`PR #142`: https://github.com/desihub/fastspecfit/pull/142 2.1.2 (2023-04-01) ------------------ diff --git a/doc/fastspec.rst b/doc/fastspec.rst index db6f8fb0..8ea56942 100644 --- a/doc/fastspec.rst +++ b/doc/fastspec.rst @@ -163,8 +163,9 @@ Name Type Units Descript FHBETA_CONT float32 1e-17 erg / (Angstrom cm2 s) Continuum flux at 4862.683 A in the rest-frame. FOIII_5007_CONT float32 1e-17 erg / (Angstrom cm2 s) Continuum flux at 5008.239 A in the rest-frame. FHALPHA_CONT float32 1e-17 erg / (Angstrom cm2 s) Continuum flux at 6564.613 A in the rest-frame. - RCHI2_LINE float32 Reduced chi-squared of an emission-line model fit. - DELTA_LINERCHI2 float32 Difference in the reduced chi-squared values between an emission-line model with narrow lines only and a model with both broad and narrow lines. + RCHI2_LINE float32 Reduced chi-squared of the emission-line model fit. + DELTA_LINECHI2 float32 Chi-squared difference between an emission-line model without and with broad lines. + DELTA_LINENDOF int32 Difference in the degrees of freedom between an emission-line model without and with broad lines. APERCORR float32 Median aperture correction factor. APERCORR_G float32 Median aperture correction factor measured in the g-band. APERCORR_R float32 Median aperture correction factor measured in the r-band. diff --git a/py/fastspecfit/emlines.py b/py/fastspecfit/emlines.py index 56104eef..9fab31e9 100644 --- a/py/fastspecfit/emlines.py +++ b/py/fastspecfit/emlines.py @@ -37,6 +37,7 @@ def build_emline_model(dlog10wave, redshift, lineamps, linevshifts, linesigmas, from fastspecfit.util import trapz_rebin #log10model = np.zeros_like(log10wave) # [erg/s/cm2/A, observed frame] + log10wave = [] # initialize # Cut to lines with non-zero amplitudes. #I = linesigmas > 0 @@ -85,14 +86,21 @@ def build_emline_model(dlog10wave, redshift, lineamps, linevshifts, linesigmas, emlinemodel = [] if camerapix is None: for icam, specwave in enumerate(emlinewave): - _emlinemodel = trapz_rebin(10**log10wave, log10model, specwave) - _emlinemodel = resolution_matrix[icam].dot(_emlinemodel) + if len(log10wave) > 0: + _emlinemodel = trapz_rebin(10**log10wave, log10model, specwave) + _emlinemodel = resolution_matrix[icam].dot(_emlinemodel) + else: + _emlinemodel = specwave * 0 emlinemodel.append(_emlinemodel) return emlinemodel else: for icam, campix in enumerate(camerapix): - _emlinemodel = trapz_rebin(10**log10wave, log10model, emlinewave[campix[0]:campix[1]]) - _emlinemodel = resolution_matrix[icam].dot(_emlinemodel) + # can be zero-length if all lines are dropped + if len(log10wave) > 0: + _emlinemodel = trapz_rebin(10**log10wave, log10model, emlinewave[campix[0]:campix[1]]) + _emlinemodel = resolution_matrix[icam].dot(_emlinemodel) + else: + _emlinemodel = emlinewave[campix[0]:campix[1]] * 0 emlinemodel.append(_emlinemodel) return np.hstack(emlinemodel) @@ -186,7 +194,7 @@ def __init__(self, minspecwave=3500.0, maxspecwave=9900.0, targetid=None): # default line-sigma for computing upper limits self.limitsigma_narrow = 75.0 self.limitsigma_broad = 1200.0 - self.wavepad = 5.0 # Angstrom + self.wavepad = 2.5 # Angstrom # Establish the names of the parameters and doublets here, at # initialization, because we use them when instantiating the best-fit @@ -212,13 +220,18 @@ def __init__(self, minspecwave=3500.0, maxspecwave=9900.0, targetid=None): self.amp_param_bool = np.array(['_amp' in pp for pp in self.param_names]) self.amp_balmer_bool = np.array(['_amp' in pp and 'hei_' not in pp and 'heii_' not in pp for pp in self.param_names]) # just Balmer lines self.sigma_param_bool = np.array(['_sigma' in pp for pp in self.param_names]) + self.vshift_param_bool = np.array(['_vshift' in pp for pp in self.param_names]) + + #self.sigma_param_bool2 = np.zeros(len(self.param_names), bool) + #self.vshift_param_bool2 = np.zeros(len(self.param_names), bool) + #for pp in self.param_names[self.amp_param_bool]: + # self.sigma_param_bool2[np.where(pp.replace('_amp', '_sigma') == self.param_names)[0]] = True + # self.vshift_param_bool2[np.where(pp.replace('_amp', '_vshift') == self.param_names)[0]] = True self.doubletindx = np.hstack([np.where(self.param_names == doublet)[0] for doublet in doublet_names]) self.doubletpair = np.hstack([np.where(self.param_names == pair)[0] for pair in doublet_pairs]) - self.delta_linerchi2_cut = 0.0 self.minsigma_balmer_broad = 250.0 # minimum broad-line sigma [km/s] - self.minsnr_balmer_broad = 3.0 def build_linemodels(self, redshift, wavelims=[3000, 10000], verbose=False, strict_finalmodel=True): """Build all the multi-parameter emission-line models we will use. @@ -751,7 +764,7 @@ def _linemodel_to_parameters(linemodel, fit_linetable): return parameters, parameter_extras @staticmethod - def populate_linemodel(linemodel, initial_guesses, param_bounds, log, broadbalmer_snrmin=3.0): + def populate_linemodel(linemodel, initial_guesses, param_bounds, log): """Population an input linemodel with initial guesses and parameter bounds, taking into account fixed parameters. @@ -770,6 +783,7 @@ def populate_linemodel(linemodel, initial_guesses, param_bounds, log, broadbalme civarkey = linemodel['linename'][iparam]+'_civar' linemodel['civar'][iparam] = initial_guesses[civarkey] + #broadbalmer_snrmin = 3. #linemodel['bounds'][iparam][0] = broadbalmer_snrmin * initial_guesses[noisekey] #if linemodel['initial'][iparam] < linemodel['bounds'][iparam][0]: # linemodel['initial'][iparam] = linemodel['bounds'][iparam][0] @@ -870,9 +884,25 @@ def optimize(self, linemodel, emlinewave, emlineflux, weights, redshift, # line-amplitude is dropped, too (see MgII 2796 on # sv1-bright-17680-39627622543528153). drop2 = np.zeros(len(parameters), bool) - if False: - drop2[Ifree] = parameters[Ifree] == linemodel['value'][Ifree] - drop2 *= notfixed + + amp_param_bool = self.amp_param_bool[Ifree] + I = np.where(parameters[Ifree][amp_param_bool] == 0)[0] + if len(I) > 0: + _Ifree = np.zeros(len(parameters), bool) + _Ifree[Ifree] = True + for pp in linemodel[Ifree][amp_param_bool][I]['param_name']: + J = np.where(_Ifree * (linemodel['param_name'] == pp.replace('_amp', '_sigma')))[0] + if len(J) > 0: + drop2[J] = True + K = np.where(_Ifree * (linemodel['param_name'] == pp.replace('_amp', '_vshift')))[0] + if len(K) > 0: + drop2[K] = True + #print(pp, J, K, np.sum(drop2)) + #pdb.set_trace() + + #drop2[self.amp_param_bool] = parameters[self.amp_param_bool] == 0.0 + #drop2[Ifree] = parameters[Ifree] == linemodel['value'][Ifree] + #drop2 *= notfixed sigmadropped = np.where(self.sigma_param_bool * drop2)[0] if len(sigmadropped) > 0: @@ -885,6 +915,17 @@ def optimize(self, linemodel, emlinewave, emlineflux, weights, redshift, drop2[linemodel['param_name'] == '{}_amp'.format(tiedline)] = True drop2[linemodel['param_name'] == '{}_amp'.format(dropline)] = True + vshiftdropped = np.where(self.vshift_param_bool * drop2)[0] + if len(vshiftdropped) > 0: + for lineindx, dropline in zip(vshiftdropped, linemodel[vshiftdropped]['linename']): + # Check whether lines are tied to this line. If so, find the + # corresponding amplitude and drop that, too. + T = linemodel['tiedtoparam'] == lineindx + if np.any(T): + for tiedline in set(linemodel['linename'][T]): + drop2[linemodel['param_name'] == '{}_amp'.format(tiedline)] = True + drop2[linemodel['param_name'] == '{}_amp'.format(dropline)] = True + # It's OK for parameters to be *at* their bounds. drop3 = np.zeros(len(parameters), bool) drop3[Ifree] = np.logical_or(parameters[Ifree] < linemodel['bounds'][Ifree, 0], @@ -892,11 +933,10 @@ def optimize(self, linemodel, emlinewave, emlineflux, weights, redshift, drop3 *= notfixed log.debug('Dropping {} negative-amplitude lines.'.format(np.sum(drop1))) # linewidth can't be negative - log.debug('Dropping {} parameters which were not optimized.'.format(np.sum(drop2))) + log.debug('Dropping {} sigma,vshift parameters of zero-amplitude lines.'.format(np.sum(drop2))) log.debug('Dropping {} parameters which are out-of-bounds.'.format(np.sum(drop3))) Idrop = np.where(np.logical_or.reduce((drop1, drop2, drop3)))[0] - - ## If we are fitting the broad Balmer lines, do some additional sanity checks. + if len(Idrop) > 0: log.debug(' Dropping {} unique parameters.'.format(len(Idrop))) parameters[Idrop] = 0.0 @@ -995,7 +1035,7 @@ def chi2(linemodel, emlinewave, emlineflux, emlineivar, emlinemodel, chi2 = 0.0 if return_dof: - return chi2, dof + return chi2, dof, nfree else: return chi2 @@ -1081,7 +1121,17 @@ def populate_emtable(self, result, finalfit, finalmodel, emlinewave, emlineflux, result[param['param_name'].upper()] = val dpixwave = emlinewave[1]-emlinewave[0] # pixel size [Angstrom] - + + # zero out all out-of-range lines + for oneline in self.fit_linetable[~self.fit_linetable['inrange']]: + linename = oneline['name'].upper() + #print(linename, result['{}_AMP'.format(linename)], result['{}_MODELAMP'.format(linename)], + # result['{}_SIGMA'.format(linename)], result['{}_VSHIFT'.format(linename)]) + result['{}_AMP'.format(linename)] = 0.0 + result['{}_MODELAMP'.format(linename)] = 0.0 + result['{}_VSHIFT'.format(linename)] = 0.0 + result['{}_SIGMA'.format(linename)] = 0.0 + # get continuum fluxes, EWs, and upper limits narrow_sigmas, broad_sigmas, uv_sigmas = [], [], [] narrow_redshifts, broad_redshifts, uv_redshifts = [], [], [] @@ -1104,7 +1154,7 @@ def populate_emtable(self, result, finalfit, finalmodel, emlinewave, emlineflux, linesigma_ang = linesigma * linezwave / C_LIGHT # [observed-frame Angstrom] - # require at least 3 pixels + # require at least 2 pixels if linesigma_ang < 2 * dpixwave: linesigma_ang_window = 2 * dpixwave else: @@ -1164,7 +1214,7 @@ def populate_emtable(self, result, finalfit, finalmodel, emlinewave, emlineflux, # np.array([result['{}_VSHIFT'.format(linename)]]), np.array([result['{}_SIGMA'.format(linename)]]), # np.array([oneline['restwave']]), emlinewave, resolution_matrix, camerapix) # - #weight = np.sum(lineprofile[lineindx]) + #weight = np.sum(lineprofile[lineindx])#**2 #if weight == 0.0: # errmsg = 'Line-profile should never sum to zero!' # log.critical(errmsg) @@ -1249,6 +1299,8 @@ def populate_emtable(self, result, finalfit, finalmodel, emlinewave, emlineflux, print() #log.debug(' ') + #if linename == 'HALPHA_BROAD': + # pdb.set_trace() ## simple QA #if linename == 'OIII_5007': # import matplotlib.pyplot as plt @@ -1584,8 +1636,6 @@ def qa_fastspec(self, data, fastspec, metadata, coadd_type='healpix', leg_uv['mgii_doublet'] = 'MgII $\lambda2796/\lambda2803={:.3f}$'.format(fastspec['MGII_DOUBLET_RATIO']) leg_broad['linerchi2'] = '$\\chi^{{2}}_{{\\nu,\\rm line}}$={:.2f}'.format(fastspec['RCHI2_LINE']) - #leg_broad['deltarchi2'] = '$\\chi^{{2}}_{{\\nu,\\rm narrow}}-\\chi^{{2}}_{{\\nu,\\rm narrow+broad}}={:.3f}$'.format(fastspec['DELTA_LINERCHI2']) - #if fastspec['DELTA_LINERCHI2'] != 0: leg_broad['deltarchi2'] = '$\\Delta\\chi^{{2}}_{{\\nu,\\rm nobroad}}={:.2f}$'.format(fastspec['DELTA_LINERCHI2']) # choose one broad Balmer line @@ -1849,7 +1899,6 @@ def qa_fastspec(self, data, fastspec, metadata, coadd_type='healpix', if np.max(modelflux) > spec_ymax: spec_ymax = np.max(modelflux) * 1.25 #print(spec_ymin, spec_ymax) - #pdb.set_trace() #specax.fill_between(wave, flux-sigma, flux+sigma, color=col1[icam], alpha=0.2) specax.plot(wave/1e4, flux, color=col1[icam], alpha=0.8) @@ -2161,8 +2210,6 @@ def major_formatter(x, pos): if _line_ymin < line_ymin[iax]: line_ymin[iax] = _line_ymin #print(linename, line_ymin[iax], line_ymax[iax]) - #if linename == '[OII] $\lambda\lambda$3726,29': - # pdb.set_trace() xx.set_xlim(wmin/1e4, wmax/1e4) @@ -2342,7 +2389,7 @@ def emline_specfit(data, templatecache, result, continuummodel, smooth_continuum modelflux """ - from astropy.table import Column + from astropy.table import Column, vstack from fastspecfit.util import ivar2var tall = time.time() @@ -2391,8 +2438,7 @@ def emline_specfit(data, templatecache, result, continuummodel, smooth_continuum initial_guesses, param_bounds = EMFit.initial_guesses_and_bounds(data, emlinewave, emlineflux, log) for linemodel in [initial_linemodel, initial_linemodel_nobroad]: - EMFit.populate_linemodel(linemodel, initial_guesses, param_bounds, log, - broadbalmer_snrmin=EMFit.minsnr_balmer_broad) + EMFit.populate_linemodel(linemodel, initial_guesses, param_bounds, log) # Initial fit - initial_linemodel_nobroad t0 = time.time() @@ -2400,105 +2446,98 @@ def emline_specfit(data, templatecache, result, continuummodel, smooth_continuum weights, redshift, resolution_matrix, camerapix, log=log, debug=False, get_finalamp=False) initmodel = EMFit.bestfit(initfit, redshift, emlinewave, resolution_matrix, camerapix) - initchi2 = EMFit.chi2(initfit, emlinewave, emlineflux, emlineivar, initmodel) - nfree = np.sum((initfit['fixed'] == False) * (initfit['tiedtoparam'] == -1)) + initchi2, initndof, nfree = EMFit.chi2(initfit, emlinewave, emlineflux, emlineivar, initmodel, return_dof=True) log.info('Initial line-fitting with {} free parameters took {:.2f} seconds [niter={}, rchi2={:.4f}].'.format( nfree, time.time()-t0, initfit.meta['nfev'], initchi2)) - # Now try adding bround Balmer and helium lines and see if we improve - # the chi2. - - ## First, do we have enough pixels around Halpha and Hbeta to do this test? - #broadlinepix = [] - #for icam in np.arange(len(data['cameras'])): - # pixoffset = int(np.sum(data['npixpercamera'][:icam])) - # for linename, linepix in zip(data['linename'][icam], data['linepix'][icam]): - # if linename == 'halpha_broad' or linename == 'hbeta_broad' or linename == 'hgamma_broad' or linename == 'hdelta_broad': - # broadlinepix.append(linepix + pixoffset) - # #print(data['wave'][icam][linepix]) - - ## Check the velocity width in the narrow lines. - #minsigma = 250.0 # [km/s] - #OIII = initfit['param_name'] == 'oiii_5007_sigma' - #if initfit['value'][OIII][0] != initfit['initial'][OIII][0]: - # if initfit['value'][OIII][0] > minsigma: - # candidate_broadline = True - # log.info('Candidate broad-line: [OIII] sigma {:.2f} > {:.0f} km/s'.format(initfit['value'][OIII][0], minsigma)) - # else: - # candidate_broadline = False - #else: - # candidate_broadline = True - - # Require minimum XX pixels and a minimum - if broadlinefit:# and candidate_broadline: - #if broadlinefit:# and len(broadlinepix) > 0 and len(np.hstack(broadlinepix)) > 10: - - ## Copy over the doublet ratios but none of the other initial values - ## so that we're not driven to the same local minimum. - #doubletindx = initial_linemodel['doubletpair'] != -1 - #doubletpair = initial_linemodel['doubletpair'][doubletindx].data - #initial_linemodel['value'][doubletindx] = initfit[doubletindx]['value'] - #initial_linemodel['value'][initial_linemodel['doubletpair'][doubletindx]] *= initial_linemodel['value'][doubletindx] - - t0 = time.time() - broadfit = EMFit.optimize(initial_linemodel, emlinewave, emlineflux, weights, - redshift, resolution_matrix, camerapix, log=log, - debug=False, get_finalamp=True) - broadmodel = EMFit.bestfit(broadfit, redshift, emlinewave, resolution_matrix, camerapix) - broadchi2 = EMFit.chi2(broadfit, emlinewave, emlineflux, emlineivar, broadmodel) - nfree = np.sum((broadfit['fixed'] == False) * (broadfit['tiedtoparam'] == -1)) - log.info('Second (broad) line-fitting with {} free parameters took {:.2f} seconds [niter={}, rchi2={:.4f}].'.format( - nfree, time.time()-t0, broadfit.meta['nfev'], broadchi2)) - linechi2_broad, linechi2_init = broadchi2, initchi2 - - log.info('Chi2 with broad lines = {:.5f} and without broad lines = {:.5f} [chi2_narrow-chi2_broad={:.5f}]'.format( - linechi2_broad, linechi2_init, linechi2_init - linechi2_broad)) - - # Choose narrow-only model if: - # --chi2_broad > chi2_narrow; - # --broad_sigma < narrow_sigma; - # --broad_sigma < 250; - # --the two reddest broad Balmer lines are both dropped. - Bbroad = broadfit['isbalmer'] * broadfit['isbroad'] * (broadfit['fixed'] == False) * EMFit.amp_balmer_bool - Habroad = broadfit['param_name'] == 'halpha_broad_sigma' - - dchi2fail = (linechi2_init - linechi2_broad) < EMFit.delta_linerchi2_cut - sigdrop1 = (broadfit[Habroad]['value'] <= broadfit[broadfit['param_name'] == 'halpha_sigma']['value'])[0] - sigdrop2 = broadfit[Habroad]['value'][0] < EMFit.minsigma_balmer_broad - - ampsnr = broadfit[Bbroad]['obsvalue'].data * np.sqrt(broadfit[Bbroad]['civar'].data) - #ampdrop = np.any(ampsnr[-1:] < EMFit.minsnr_balmer_broad) - ampdrop = np.any(ampsnr[-2:] < EMFit.minsnr_balmer_broad) - - #W = (initfit['fixed'] == False) * (initfit['tiedtoparam']==-1) - #W2 = (broadfit['fixed'] == False) * (broadfit['tiedtoparam']==-1) - - if dchi2fail or ampdrop or sigdrop1 or sigdrop2: - if dchi2fail: - log.info('Dropping broad-line model: delta-rchi2 {:.3f}<{:.3f}.'.format(linechi2_init - linechi2_broad, EMFit.delta_linerchi2_cut)) - elif ampdrop: - log.info('Dropping broad-line model: S/N in either of the two reddest broad lines < {:.1f}.'.format(EMFit.minsnr_balmer_broad)) - #if alldrop: - # log.info('Dropping broad-line model: all broad lines dropped.') - elif sigdrop1: - log.info('Dropping broad-line model: Halpha_broad_sigma {:.2f} km/s < Halpha_narrow_sigma {:.2f} km/s.'.format( - broadfit[Habroad]['value'][0], - broadfit[broadfit['param_name'] == 'halpha_sigma']['value'][0])) - elif sigdrop2: - log.info('Dropping broad-line model: Halpha_broad_sigma {:.2f} km/s < {:.0f} km/s.'.format( - broadfit[Habroad]['value'][0], EMFit.minsigma_balmer_broad)) + # Now try adding broad Balmer and helium lines and see if we improve the + # chi2. + if broadlinefit: + # Gather the pixels around the broad Balmer lines and the corresponding + # linemodel table. + balmer_pix, balmer_linemodel, balmer_linemodel_nobroad = [], [], [] + for icam in np.arange(len(data['cameras'])): + pixoffset = int(np.sum(data['npixpercamera'][:icam])) + if len(data['linename'][icam]) > 0: + I = (initial_linemodel['isbalmer'] * initial_linemodel['isbroad'] * + np.isin(initial_linemodel['linename'], data['linename'][icam])) + _balmer_linemodel = initial_linemodel[I] + _balmer_linemodel_nobroad = initial_linemodel_nobroad[I] + balmer_linemodel.append(_balmer_linemodel) + balmer_linemodel_nobroad.append(_balmer_linemodel_nobroad) + if len(_balmer_linemodel) > 0: # use balmer_linemodel not balmer_linemodel_nobroad + I = np.where(np.isin(data['linename'][icam], _balmer_linemodel['linename']))[0] + for ii in I: + #print(data['linename'][icam][ii]) + balmer_pix.append(data['linepix'][icam][ii] + pixoffset) + + if len(balmer_pix) > 0: + t0 = time.time() + broadfit = EMFit.optimize(initial_linemodel, emlinewave, emlineflux, weights, + redshift, resolution_matrix, camerapix, log=log, + debug=False, get_finalamp=True) + broadmodel = EMFit.bestfit(broadfit, redshift, emlinewave, resolution_matrix, camerapix) + broadchi2, broadndof, nfree = EMFit.chi2(broadfit, emlinewave, emlineflux, emlineivar, broadmodel, return_dof=True) + log.info('Second (broad) line-fitting with {} free parameters took {:.2f} seconds [niter={}, rchi2={:.4f}].'.format( + nfree, time.time()-t0, broadfit.meta['nfev'], broadchi2)) + + # compute delta-chi2 around just the Balmer lines + balmer_pix = np.hstack(balmer_pix) + balmer_linemodel = vstack(balmer_linemodel) + + balmer_nfree = (np.count_nonzero((balmer_linemodel['fixed'] == False) * + (balmer_linemodel['tiedtoparam'] == -1))) + balmer_ndof = np.count_nonzero(emlineivar[balmer_pix] > 0) - balmer_nfree + + balmer_linemodel_nobroad = vstack(balmer_linemodel_nobroad) + balmer_nfree_nobroad = (np.count_nonzero((balmer_linemodel_nobroad['fixed'] == False) * + (balmer_linemodel_nobroad['tiedtoparam'] == -1))) + balmer_ndof_nobroad = np.count_nonzero(emlineivar[balmer_pix] > 0) - balmer_nfree_nobroad + + linechi2_balmer = np.sum(emlineivar[balmer_pix] * (emlineflux[balmer_pix] - broadmodel[balmer_pix])**2) + linechi2_balmer_nobroad = np.sum(emlineivar[balmer_pix] * (emlineflux[balmer_pix] - initmodel[balmer_pix])**2) + delta_linechi2_balmer = linechi2_balmer_nobroad - linechi2_balmer + delta_linendof_balmer = balmer_ndof_nobroad - balmer_ndof + + # Choose broad-line model only if: + # --delta-chi2 > delta-ndof + # --broad_sigma < narrow_sigma + # --broad_sigma < 250 + + dchi2test = delta_linechi2_balmer > delta_linendof_balmer + Hanarrow = broadfit['param_name'] == 'halpha_sigma' # Balmer lines are tied to H-alpha even if out of range + Habroad = broadfit['param_name'] == 'halpha_broad_sigma' + sigtest1 = broadfit[Habroad]['value'][0] > EMFit.minsigma_balmer_broad + sigtest2 = (broadfit[Habroad]['value'] > broadfit[Hanarrow]['value'])[0] + + if dchi2test and sigtest1 and sigtest2: + log.info('Adopting broad-line model: delta-chi2={:.1f} > delta-ndof={:.0f} (sigma_broad={:.1f} km/s, sigma_narrow={:.1f} km/s)'.format( + delta_linechi2_balmer, delta_linendof_balmer, broadfit[Habroad]['value'][0], broadfit[Hanarrow]['value'][0])) + bestfit = broadfit + use_linemodel_broad = True + else: + if dchi2test == False: + log.info('Dropping broad-line model: delta-chi2={:.1f} < delta-ndof={:.0f}'.format( + delta_linechi2_balmer, delta_linendof_balmer)) + elif sigtest1 == False: + log.info('Dropping broad-line model: Halpha_broad_sigma {:.1f} km/s < {:.0f} km/s (delta-chi2={:.1f}, delta-ndof={:.0f}).'.format( + broadfit[Habroad]['value'][0], EMFit.minsigma_balmer_broad, delta_linechi2_balmer, delta_linendof_balmer)) + elif sigtest2 == False: + log.info('Dropping broad-line model: Halpha_broad_sigma {:.1f} km/s < Halpha_narrow_sigma {:.1f} km/s (delta-chi2={:.1f}, delta-ndof={:.0f}).'.format( + broadfit[Habroad]['value'][0], broadfit[Hanarrow]['value'][0], delta_linechi2_balmer, delta_linendof_balmer)) + bestfit = initfit + use_linemodel_broad = False + else: + log.info('Insufficient Balmer lines to test the broad-line model.') bestfit = initfit + delta_linechi2_balmer, delta_linendof_balmer = 0, np.int32(0) use_linemodel_broad = False - else: - log.info('Adopting broad-line model: delta-rchi2={:.3f}>{:.3f}'.format(linechi2_init - linechi2_broad, EMFit.delta_linerchi2_cut)) - bestfit = broadfit - use_linemodel_broad = True else: log.info('Skipping broad-line fitting (broadlinefit=False).') bestfit = initfit - linechi2_broad, linechi2_init = 1e6, initchi2 + delta_linechi2_balmer, delta_linendof_balmer = 0, np.int32(0) use_linemodel_broad = False - + # Finally, one more fitting loop with all the line-constraints relaxed # but starting from the previous best-fitting values. if use_linemodel_broad: @@ -2581,8 +2620,7 @@ def emline_specfit(data, templatecache, result, continuummodel, smooth_continuum redshift, resolution_matrix, camerapix, log=log, debug=False, get_finalamp=True) finalmodel = EMFit.bestfit(finalfit, redshift, emlinewave, resolution_matrix, camerapix) - finalchi2 = EMFit.chi2(finalfit, emlinewave, emlineflux, emlineivar, finalmodel) - nfree = np.sum((finalfit['fixed'] == False) * (finalfit['tiedtoparam'] == -1)) + finalchi2, finalndof, nfree = EMFit.chi2(finalfit, emlinewave, emlineflux, emlineivar, finalmodel, return_dof=True) log.info('Final line-fitting with {} free parameters took {:.2f} seconds [niter={}, rchi2={:.4f}].'.format( nfree, time.time()-t0, finalfit.meta['nfev'], finalchi2)) @@ -2607,7 +2645,9 @@ def emline_specfit(data, templatecache, result, continuummodel, smooth_continuum emmodel = np.hstack(EMFit.emlinemodel_bestfit(data['wave'], data['res'], result, redshift=redshift)) result['RCHI2_LINE'] = finalchi2 - result['DELTA_LINERCHI2'] = linechi2_init - linechi2_broad + #result['NDOF_LINE'] = finalndof + result['DELTA_LINECHI2'] = delta_linechi2_balmer # linechi2_init - linechi2_broad + result['DELTA_LINENDOF'] = delta_linendof_balmer # linendof_init - linendof_broad # full-fit reduced chi2 rchi2 = np.sum(oemlineivar * (specflux - (continuummodelflux + smoothcontinuummodelflux + emmodel))**2) diff --git a/py/fastspecfit/fastspecfit.py b/py/fastspecfit/fastspecfit.py index d77d076e..d4baaaac 100644 --- a/py/fastspecfit/fastspecfit.py +++ b/py/fastspecfit/fastspecfit.py @@ -523,7 +523,9 @@ def major_formatter(x, pos): leg_broad['linerchi2'] = '$\\chi^{{2}}_{{\\nu,\\rm line}}$={:.2f}'.format(fastspec['RCHI2_LINE']) #leg_broad['deltarchi2'] = '$\\chi^{{2}}_{{\\nu,\\rm narrow}}-\\chi^{{2}}_{{\\nu,\\rm narrow+broad}}={:.3f}$'.format(fastspec['DELTA_LINERCHI2']) #if fastspec['DELTA_LINERCHI2'] != 0: - leg_broad['deltarchi2'] = '$\\Delta\\chi^{{2}}_{{\\nu,\\rm nobroad}}={:.2f}$'.format(fastspec['DELTA_LINERCHI2']) + #leg_broad['deltachi2'] = '$\\Delta\\chi^{{2}}_{{\\nu,\\rm nobroad}}={:.2f}$'.format(fastspec['DELTA_LINECHI2']) + leg_broad['deltachi2'] = '$\\Delta\\chi^{{2}}_{{\\rm nobroad}}={:.2f}$'.format(fastspec['DELTA_LINECHI2']) + leg_broad['deltandof'] = '$\\Delta\\nu_{{\\rm nobroad}}={:.0f}$'.format(fastspec['DELTA_LINENDOF']) # choose one broad Balmer line if fastspec['HALPHA_BROAD_AMP']*np.sqrt(fastspec['HALPHA_BROAD_AMP_IVAR']) > snrcut: diff --git a/py/fastspecfit/io.py b/py/fastspecfit/io.py index ab52a18a..4d07b0d9 100644 --- a/py/fastspecfit/io.py +++ b/py/fastspecfit/io.py @@ -19,10 +19,10 @@ log = get_logger() # Default environment variables. -DESI_ROOT_NERSC = '/global/cfs/cdirs/desi' -DUST_DIR_NERSC = '/global/cfs/cdirs/cosmo/data/dust/v0_1' -DR9_DIR_NERSC = '/global/cfs/cdirs/desi/external/legacysurvey/dr9' -FTEMPLATES_DIR_NERSC = '/global/cfs/cdirs/desi/science/gqp/templates/fastspecfit' +DESI_ROOT_NERSC = '/dvs_ro/cfs/cdirs/desi' +DUST_DIR_NERSC = '/dvs_ro/cfs/cdirs/cosmo/data/dust/v0_1' +DR9_DIR_NERSC = '/dvs_ro/cfs/cdirs/desi/external/legacysurvey/dr9' +FTEMPLATES_DIR_NERSC = '/dvs_ro/cfs/cdirs/desi/science/gqp/templates/fastspecfit' # list of all possible targeting bit columns TARGETINGBITS = { @@ -1634,8 +1634,10 @@ def init_fastspec_output(input_meta, specprod, templates=None, ncoeff=None, # Add chi2 metrics #out.add_column(Column(name='DOF', length=nobj, dtype='i8')) # full-spectrum dof out.add_column(Column(name='RCHI2_LINE', length=nobj, dtype='f4')) # reduced chi2 with broad line-emission + #out.add_column(Column(name='NDOF_LINE', length=nobj, dtype='i8')) # number of degrees of freedom corresponding to rchi2_line #out.add_column(Column(name='DOF_BROAD', length=nobj, dtype='i8')) - out.add_column(Column(name='DELTA_LINERCHI2', length=nobj, dtype='f4')) # delta-reduced chi2 with and without broad line-emission + out.add_column(Column(name='DELTA_LINECHI2', length=nobj, dtype='f4')) # delta-reduced chi2 with and without broad line-emission + out.add_column(Column(name='DELTA_LINENDOF', length=nobj, dtype=np.int32)) # aperture corrections out.add_column(Column(name='APERCORR', length=nobj, dtype='f4')) # median aperture correction