diff --git a/py/fastspecfit/continuum.py b/py/fastspecfit/continuum.py index 14bc405b..419d32d1 100644 --- a/py/fastspecfit/continuum.py +++ b/py/fastspecfit/continuum.py @@ -1422,6 +1422,7 @@ def init_spec_output(self, nobj=1): nssp_coeff = len(self.sspinfo) out = Table() + out.add_column(Column(name='APCORR', length=nobj, dtype='f4')) # aperture correction out.add_column(Column(name='CONTINUUM_Z', length=nobj, dtype='f8')) # redshift out.add_column(Column(name='CONTINUUM_COEFF', length=nobj, shape=(nssp_coeff,), dtype='f8')) out.add_column(Column(name='CONTINUUM_RCHI2', length=nobj, dtype='f4')) # reduced chi2 @@ -2126,8 +2127,8 @@ def continuum_specfit(self, data): bestfit = bestsspflux.dot(coeff) meanage = self.get_meanage(coeff) - flam_ivar = np.hstack(data['ivar']) - dn4000, dn4000_ivar = self.get_dn4000(specwave, specflux, flam_ivar=flam_ivar, # specivar is line-masked! + flam_ivar = np.hstack(data['ivar']) # specivar is line-masked! + dn4000, dn4000_ivar = self.get_dn4000(specwave, specflux, flam_ivar=flam_ivar, redshift=redshift, rest=False) dn4000_model, _ = self.get_dn4000(specwave, bestfit, redshift=redshift, rest=False) @@ -2284,6 +2285,7 @@ def continuum_joint_specfit(self, data): # smooth continuum correction), which will guarantee that the aperture # correction is always positive. apcorr = data['phot']['nanomaggies'][1] / data['synthphot']['nanomaggies'][1] + data['apcorr'] = apcorr # Prepare the reddened and unreddened SSP templates by redshifting and # normalizing. Note that we ignore templates which are older than the @@ -2442,24 +2444,28 @@ def continuum_joint_specfit(self, data): self.log.info('Insufficient S/N to compute vdisp; adopting vdisp={:.2f} km/s'.format(self.vdisp_nominal)) vdispbest, vdispivar = self.vdisp_nominal, 0.0 - pdb.set_trace() - # Get the final set of coefficients and chi2 at the best-fitting # reddening and velocity dispersion. - bestsspflux, _ = self.SSP2data(self.sspflux[:, agekeep], self.sspwave, redshift=redshift, - specwave=data['wave'], specres=data['res'], - AV=AVbest, vdisp=vdispbest, cameras=data['cameras'], - south=data['photsys'] == 'S', synthphot=False) + # dev - synthesize photometry + bestsspflux, bestphot = self.SSP2data(self.sspflux[:, agekeep], self.sspwave, redshift=redshift, + specwave=data['wave'], specres=data['res'], + AV=AVbest, vdisp=vdispbest, cameras=data['cameras'], + south=data['photsys'] == 'S', synthphot=True) bestsspflux = np.concatenate(bestsspflux, axis=0) - coeff, chi2min = self._call_nnls(bestsspflux, specflux, specivar) - chi2min /= np.sum(specivar > 0) # dof??? + bestflam = bestphot['flam'].data*self.massnorm*self.fluxnorm + + # dev - prepend the photometry + coeff, chi2min = self._call_nnls(np.vstack((bestflam, bestsspflux)), + np.hstack((objflam, specflux*apcorr)), + np.hstack((objflamivar, specivar/apcorr**2))) + chi2min /= (np.sum(objflamivar > 0) + np.sum(specivar > 0)) # dof??? # Get the light-weighted age and DN(4000). bestfit = bestsspflux.dot(coeff) meanage = self.get_meanage(coeff) - flam_ivar = np.hstack(data['ivar']) - dn4000, dn4000_ivar = self.get_dn4000(specwave, specflux, flam_ivar=flam_ivar, # specivar is line-masked! + flam_ivar = np.hstack(data['ivar']) # specivar is line-masked! + dn4000, dn4000_ivar = self.get_dn4000(specwave, specflux, flam_ivar=flam_ivar, redshift=redshift, rest=False) dn4000_model, _ = self.get_dn4000(specwave, bestfit, redshift=redshift, rest=False) @@ -2511,7 +2517,7 @@ def continuum_joint_specfit(self, data): if np.all(coeff == 0): _smooth_continuum = np.zeros_like(bestfit) else: - _smooth_continuum, _ = self.smooth_continuum(specwave, specflux - bestfit, specivar, + _smooth_continuum, _ = self.smooth_continuum(specwave, apcorr*specflux - bestfit, specivar/apcorr**2, redshift, linemask=linemask, png=png) # Unpack the continuum into individual cameras. @@ -2522,24 +2528,25 @@ def continuum_joint_specfit(self, data): ## Like above, but with per-camera smoothing. #smooth_continuum = self.smooth_residuals( - # continuummodel, data['wave'], data['flux'], + # continuummodel, data['wave'], apcorr*data['flux'], # data['ivar'], data['linemask'], percamera=False) if False: import matplotlib.pyplot as plt fig, ax = plt.subplots(2, 1) for icam in np.arange(len(data['cameras'])): # iterate over cameras - resid = data['flux'][icam]-continuummodel[icam] + resid = apcorr*data['flux'][icam]-continuummodel[icam] ax[0].plot(data['wave'][icam], resid) ax[1].plot(data['wave'][icam], resid-smooth_continuum[icam]) for icam in np.arange(len(data['cameras'])): # iterate over cameras - resid = data['flux'][icam]-continuummodel[icam] + resid = apcorr*data['flux'][icam]-continuummodel[icam] pix_emlines = np.logical_not(data['linemask'][icam]) # affected by line = True ax[0].scatter(data['wave'][icam][pix_emlines], resid[pix_emlines], s=30, color='red') ax[0].plot(data['wave'][icam], smooth_continuum[icam], color='k', alpha=0.7, lw=2) plt.savefig('junk.png') # Pack it in and return. + result['APCORR'] = apcorr result['CONTINUUM_COEFF'][0][0:nage] = coeff result['CONTINUUM_RCHI2'][0] = chi2min result['CONTINUUM_AV'][0] = AVbest # * u.mag diff --git a/py/fastspecfit/emlines.py b/py/fastspecfit/emlines.py index 24f22188..68137c8d 100644 --- a/py/fastspecfit/emlines.py +++ b/py/fastspecfit/emlines.py @@ -981,8 +981,8 @@ def fit(self, data, continuummodel, smooth_continuum, synthphot=True, # best-fitting model (per-camera) below. redshift = data['zredrock'] emlinewave = np.hstack(data['wave']) - oemlineivar = np.hstack(data['ivar']) - specflux = np.hstack(data['flux']) + oemlineivar = np.hstack(data['ivar'])/data['apcorr']**2 + specflux = np.hstack(data['flux'])*data['apcorr'] resolution_matrix = data['res'] camerapix = data['camerapix'] @@ -1593,7 +1593,9 @@ def qa_fastspec(self, data, fastspec, metadata, coadd_type='healpix', else: pass + apcorr = fastspec['APCORR'] redshift = fastspec['CONTINUUM_Z'] + stackwave = np.hstack(data['wave']) # rebuild the best-fitting spectroscopic and photometric models @@ -1604,14 +1606,23 @@ def qa_fastspec(self, data, fastspec, metadata, coadd_type='healpix', vdisp=fastspec['CONTINUUM_VDISP'], coeff=fastspec['CONTINUUM_COEFF'], synthphot=False) - - residuals = [data['flux'][icam] - continuum[icam] for icam in np.arange(len(data['cameras']))] + print('This bit is new!') + continuum_phot, synthmodelphot = self.SSP2data( + self.sspflux, self.sspwave, redshift=redshift, + synthphot=True, AV=fastspec['CONTINUUM_AV'], #test=True, + vdisp=fastspec['CONTINUUM_VDISP'], + coeff=fastspec['CONTINUUM_COEFF'] * self.massnorm) + continuum_wave_phot = self.sspwave * (1 + redshift) + + residuals = [apcorr*data['flux'][icam] - continuum[icam] for icam in np.arange(len(data['cameras']))] + if np.all(fastspec['CONTINUUM_COEFF'] == 0): _smooth_continuum = np.zeros_like(stackwave) else: _smooth_continuum, _ = self.smooth_continuum(np.hstack(data['wave']), np.hstack(residuals), np.hstack(data['ivar']), redshift=redshift, linemask=np.hstack(data['linemask'])) + smooth_continuum = [] for campix in data['camerapix']: smooth_continuum.append(_smooth_continuum[campix[0]:campix[1]]) @@ -1697,23 +1708,23 @@ def qa_fastspec(self, data, fastspec, metadata, coadd_type='healpix', #ymin = np.zeros(len(data['cameras'])) #ymax = np.zeros(len(data['cameras'])) for ii in np.arange(len(data['cameras'])): # iterate over cameras - sigma, good = ivar2var(data['ivar'][ii], sigma=True, allmasked_ok=True) + sigma, good = ivar2var(data['ivar'][ii]/apcorr**2, sigma=True, allmasked_ok=True) - bigax1.fill_between(data['wave'][ii], data['flux'][ii]-sigma, - data['flux'][ii]+sigma, color=col1[ii]) + bigax1.fill_between(data['wave'][ii], apcorr*data['flux'][ii]-sigma, + apcorr*data['flux'][ii]+sigma, color=col1[ii]) #bigax1.plot(data['wave'][ii], continuum_nodust[ii], alpha=0.5, color='k') #bigax1.plot(data['wave'][ii], smooth_continuum[ii], color='gray')#col3[ii])#, alpha=0.3, lw=2)#, color='k') #bigax1.plot(data['wave'][ii], continuum[ii], color=col2[ii]) bigax1.plot(data['wave'][ii], continuum[ii]+smooth_continuum[ii], color=col2[ii]) # get the robust range - filtflux = median_filter(data['flux'][ii], 51, mode='nearest') + filtflux = median_filter(apcorr*data['flux'][ii], 51, mode='nearest') #filtflux = median_filter(data['flux'][ii] - _emlinemodel[ii], 51, mode='nearest') #perc = np.percentile(filtflux[data['ivar'][ii] > 0], [5, 95]) - #sigflux = np.std(data['flux'][ii][data['ivar'][ii] > 0]) + #sigflux = np.std(apcorr*data['flux'][ii][data['ivar'][ii] > 0]) I = data['ivar'][ii] > 0 if np.sum(I) > 0: - sigflux = np.diff(np.percentile(data['flux'][ii][I], [25, 75]))[0] / 1.349 # robust + sigflux = np.diff(np.percentile(apcorr*data['flux'][ii][I], [25, 75]))[0] / 1.349 # robust else: sigflux = 0.0 #sigflux = np.std(filtflux[data['ivar'][ii] > 0]) @@ -1743,6 +1754,8 @@ def qa_fastspec(self, data, fastspec, metadata, coadd_type='healpix', bigax1.plot(stackwave, _smooth_continuum, color='gray')#col3[ii])#, alpha=0.3, lw=2)#, color='k') bigax1.plot(stackwave, np.hstack(continuum), color='k', alpha=0.1)#col3[ii])#, alpha=0.3, lw=2)#, color='k') + #bigax1.plot(continuum_wave_phot, continuum_phot / self.massnorm, color='orange', lw=4) + bigax1.text(0.03, 0.9, 'Observed Spectrum + Continuum Model', ha='left', va='center', transform=bigax1.transAxes, fontsize=30) if not self.nolegend: @@ -1770,10 +1783,10 @@ def qa_fastspec(self, data, fastspec, metadata, coadd_type='healpix', ymin, ymax = 1e6, -1e6 for ii in np.arange(len(data['cameras'])): # iterate over cameras emlinewave = data['wave'][ii] - emlineflux = data['flux'][ii] - continuum[ii] - smooth_continuum[ii] + emlineflux = apcorr*data['flux'][ii] - continuum[ii] - smooth_continuum[ii] emlinemodel = _emlinemodel[ii] - emlinesigma, good = ivar2var(data['ivar'][ii], sigma=True, allmasked_ok=True) + emlinesigma, good = ivar2var(data['ivar'][ii]/apcorr**2, sigma=True, allmasked_ok=True) emlinewave = emlinewave[good] emlineflux = emlineflux[good] emlinesigma = emlinesigma[good] @@ -1884,10 +1897,10 @@ def qa_fastspec(self, data, fastspec, metadata, coadd_type='healpix', # iterate over cameras for ii in np.arange(len(data['cameras'])): # iterate over cameras emlinewave = data['wave'][ii] - emlineflux = data['flux'][ii] - continuum[ii] - smooth_continuum[ii] + emlineflux = apcorr*data['flux'][ii] - continuum[ii] - smooth_continuum[ii] emlinemodel = _emlinemodel[ii] - emlinesigma, good = ivar2var(data['ivar'][ii], sigma=True, allmasked_ok=True) + emlinesigma, good = ivar2var(data['ivar'][ii]/apcorr**2, sigma=True, allmasked_ok=True) emlinewave = emlinewave[good] emlineflux = emlineflux[good] emlinesigma = emlinesigma[good] diff --git a/py/fastspecfit/fastspecfit.py b/py/fastspecfit/fastspecfit.py index a29a8dcf..f08b243b 100644 --- a/py/fastspecfit/fastspecfit.py +++ b/py/fastspecfit/fastspecfit.py @@ -458,7 +458,7 @@ def build_webqa(CFit, EMFit, data, fastfit, metadata, coadd_type='healpix', _emlinemodel_oneline.append(_emlinemodel_oneline1) # Grab the viewer cutout. - width = int(60 / 0.262) # =1 arcmin + width = int(40 / 0.262) # =1 arcmin height = int(width / 1.3) # 3:2 aspect ratio cutoutpng = os.path.join('/tmp', 'tmp.'+os.path.basename(pngfile))