DAP IDL Tutorial
The goal of this introductory tutorial will be to show you the basics of loading and visualizing MaNGA DAP products. It assumes you're using IDL and have the IDL Astronomy User's Library, the Coyote Graphics Library, and image_plot.pro from the ppxf package. This tutorial was last tested with IDL 8.8.0.
Loading a MAPS file
For this tutorial, we'll use the maps file for MaNGA object (plate-ifu) 7443-12703 with the MAPS-HYB10-MILESHC-MASTARSSP type map which can be downloaded here.
file={dir}+'manga-7443-12703-MAPS-HYB10-MILESHC-MASTARSSP.fits.gz'
main=mrdfits(file,0,mainhdr,/silent)
Here {dir} is the directory where you downloaded the maps file. First, we read in the main header in the 0th extension. This header has some important information like the DAPQUAL flag. Let's check to make sure that this flag is 0.
qual=sxpar(mainhdr,'DAPQUAL') print,qual ;should print 0
Making an emission line map and applying masks
Now let's make an emission line map and apply the appropriate masks. First, we need to read in the emission-line flux extension, and we choose the Gaussian-fitted flux in this example instead of the summed flux, labeled EMLINE_GFLUX. The different extensions and descriptions can be found here.
emline_flux=mrdfits(file,'EMLINE_GFLUX',fluxhdr,/silent) help,emline_flux ;see the size and type of emline_flux
We see that the {emline_flux} is a float of size [74,74,35]. This IFU has 74 by 74 spaxels and there are 35 channels for the different emission lines.
We want to read in the Hα flux. First, we need to find the channel with the Hα flux among all the channels in the EMLINE_GFLUX extension. This information is in the header of the extension.
print,fluxhdr
We see in the header that "C24 = 'Ha-6564' /Data in channel 24".
ha_flux=emline_flux[*,*,23] ;since IDL starts indices with 0 help,ha_flux ; should say HA_FLUX FLOAT =Array[74,74]
For making 2D map plots in IDL, you can use 'image_plot.pro'
x=findgen(74) y=findgen(74) image_plot,ha_flux,x,y
To change the color map used in IDL, try changing the {loadct} to 35.
loadct,35 image_plot,ha_flux,x,y ;replot with the new color table
There are many spaxels with unreliable Hα flux measurements. We should apply the mask using the provided mask extension, specifically for the Hα measurements.
emline_mask=mrdfits(file,'EMLINE_GFLUX_MASK',/silent) ha_mask=emline_mask[*,*,23] ha_good=ha_flux ha_good(where(ha_mask ne 0))=-99 ;applying the mask and setting masked spaxels to flux of -99 image_plot,ha_good,x,y
Plotting a velocity field and creating your own masks
We will now follow a similar procedure as above for the Hα fluxes, but now for a velocity field.
First we need to read in the velocity extension and Hα channel.
gas_vfield=mrdfits(file,'EMLINE_GVEL',/silent) gas_vfield_mask=mrdfits(file,'EMLINE_GVEL_MASK',/silent) ha_vfield=gas_vfield[*,*,23] ha_vfield_mask=gas_vfield_mask[*,*,23] ha_vfield_good=ha_vfield ha_vfield_good(where(ha_vfield_mask ne 0))=-999 image_plot,ha_vfield_good,x,y
Let's now make a signal-to-noise cut, so only the Hα velocities with a decent signal-to-noise are shown. For this we will use the inverse variance extension.
emline_ivar=mrdfits(file,'EMLINE_GFLUX_IVAR',/silent) ha_ivar=emline_ivar[*,*,23] ha_sn=ha_flux*sqrt(ha_ivar) ;calculating the signal-to-noise for Hα sn_cut=10. ha_vfield_sn=ha_vfield_good ha_vfield_sn(where(ha_sn lt sn_cut))=-999 ;applying the cut image_plot,ha_vfield_sn,x,y
An important note regarding gas velocity fields: The velocities of the emission lines are tied together, so for example, the velocity of the [OIII]-5007 line is the same as the Hα line, as are the uncertainties. You cannot reduce the uncertainty on the measured velocity by averaging the velocities of several lines together.
Plotting stellar velocity dispersion
In this section, we will plot the stellar velocity dispersion. First we need to read in the raw stellar velocities.
disp_raw=mrdfits(file,'STELLAR_SIGMA',/silent)
We need to correct the stellar velocity dispersion for the instrumental resolution.
disp_inst=mrdfits(file,'STELLAR_SIGMACORR',/silent)
Let’s apply the correction and plot the results (also removing masked values). The calculation below will ignore any points where the correction is larger than the measured dispersion.
disp_stars_corr=sqrt(disp_raw^2-disp_inst^2) disp_mask=mrdfits(file,'STELLAR_SIGMA_MASK',/silent) ;reading in mask disp_stars_corr(where(disp_mask ne 0 or finite(disp_stars_corr) eq 0))=-99 ;applying mask and removing NaN values image_plot,disp_stars_corr,x,y
Plotting spectral-index measurements
Next we’ll show how to plot spectral indices and correct them for velocity dispersion. The spectral indices are stored in the SPECINDEX extension of the MAPS file. We will follow a similar procedure as before with the emission line fluxes.
specindex=mrdfits(file,'SPECINDEX',spechdr,/silent) print,spechdr
Like with the emission-line header, the spectral indices header provides the name of the spectral index in each channel. Also the units for each index is given, which is important for the corrections done below. We see that Hβ is in channel 9 and the units are in Angstrom (C09='Hb' and U09='ang'). We will also get the masks and the corrections for the velocity dispersion for the Hβ spectral index.
hb=specindex[*,*,8] specindex_mask=mrdfits(file,'SPECINDEX_MASK',/silent) specindex_corr=mrdfits(file,'SPECINDEX_CORR',/silent) hb_mask=specindex_mask[*,*,8] hb_corr=specindex_corr[*,*,8]
Now we can apply the velocity dispersion correction. Since Hβ spectral index is in Angstroms, we can corrected it like this:
hb_vdc=hb*hb_corr
and apply the mask like this:
hb_vdc(where(hb_mask ne 0))=-150 ;choose -150 so the plot is more legible image_plot,hb_vdc,x,y
Identifying unique bins
The spaxels are binned in different ways depending on the measurement being made (the DAP documentation provides more information). This binning means that two spaxels can belong to the same bin, and therefore a derived quantity at those locations will be identical. The BINID extension provides information about which spaxels are in each bin. There are 5 channels providing the IDs of spaxels associated with
- 0: each binned spectrum
- 1: any binned spectrum with an attempted stellar kinematics fit.
- 2: any binned spectrum with emission-line moment measurements.
- 3: any binned spectrum with an attempted emission-line fit.
- 4: any binned spectrum with spectral-index measurements.
Any spaxel in any of these channels with a value of -1 means that the spaxel was not included in the analysis, e.g., because its S/N was too low.
For any analysis, you'll want to extract the unique spectra and/or maps values. For instance, to find the indices of the unique bins where stellar kinematics were fit:
binids=mrdfits(file,'binid') stellar_kin_bins = binids[*,*,1] unique_bin_indices = uniq(stellar_kin_bins,sort(stellar_kin_bins)) unique_bins = stellar_kin_bins[unique_bin_indices]
Since the unique bins are returned as an ordered array, we can just remove the first elements to get rid of bin IDs set to -1 (meaning they were not analyzed):
n_ind = n_elements(unique_bin_indices) unique_bin_indices = unique_bin_indices[1:n_ind-1] unique_bins = unique_bins[1:n_ind-1]
Let's now use this information to plot the position of each unique bin and color-code it by the measured stellar velocity:
xy = mrdfits(file,'BIN_LWSKYCOO') ;coordinates of each spaxel
x=xy[*,*,0]
y=xy[*,*,1]
vstars = mrdfits(file,'STELLAR_VEL') ;stellar velocity
vstars_mask = mrdfits(file,'STELLAR_VEL_MASK') ;stellar velocity mask
vmax=150
vmin=-150
colors = round((vstars-vmin)/(vmax-vmin)*255) ;scale color by velocity
colors = colors > 0
colors = colors < 255
device,decomposed=0
loadct,0,/silent
cgplot,x[unique_bin_indices],y[unique_bin_indices],psym=4,$
xtitle='x (arcsec)',ytitle='y (arcsec)',/nodata,position=[.125,.125,0.8,0.9]
loadct,70,/silent
plotsym,0,1,/fill
for i=0,n_elements(unique_bin_indices)-1 do $
if vstars_mask[unique_bin_indices[i]] ne 1 then $
cgoplot,x[unique_bin_indices[i]],y[unique_bin_indices[i]],psym=8,$
color=colors[unique_bin_indices[i]]
cgcolorbar,range=[vmin,vmax],/vertical,/right,title='stellar velocity (km/s)'
Extract a binned spectrum and its model fit
The model cube files provide detailed information about the output the binned spaxels and the model fitting. We can use these files to compare an individual bin's measured spectrum, model fit, model residuals, and so on. However, one needs to be careful when comparing binned spectra with the model fits. Specifically, there are three types of files with different binning schemes:
SPX: All spaxels with S/N ≥ 1 are analyzed.VOR10: The spectra are Voronoi binned to S/N ˜ 10. Stellar and emission line parameters are estimated from those bins.HYB10: The spectra are again Voronoi binned to S/N ˜ 10, and the stellar parameters are calculated using these Voronoi binned spectra. However, emission-line parameters are measured using the individual 0.5"x0.5" spaxels.
If you are using the HYB10 files and want to compare the best fitting model (including stellar continuum and emission lines) to the data, you need to compare the models to the individual spectra (measured in 0.5"x0.5" spaxels) from the DRP LOGCUBE files, not the binned spectra in the HYB10 files. If you are comparing the best-fitting stellar continuum models to the data, you should use the binned spectra within the HYB10 files.
Below are a few examples using both VOR10 and HYB10 files to demonstrate the proper way to compare the models and data using these two files types.
VOR10 Files
VOR10 files where comparing the models and data is simpler. First, let's read in the BIN_SNR extension from the VOR10 maps file and find the bin with the highest S/N. The file for this example can be downloaded here.
file=dir+'manga-7443-12703-MAPS-VOR10-MILESHC-MASTARSSP.fits.gz' snr=mrdfits(file,'BIN_SNR') maxsnr = max(snr,ind) xy = array_indices(snr,ind) ;convert into 2d indices
Next we'll load the model cube file for this galaxy (can be downloaded here, pull out the binned spectrum at its location, and then plot the binned spectrum, the full best fit model, the model stellar continuum, the model emission lines, and the residuals.
modelcube = dir+'manga-7443-12703-LOGCUBE-VOR10-MILESHC-MASTARSSP.fits.gz'
wave = mrdfits(modelcube,'WAVE') ;wavelength array
flux = mrdfits(modelcube,'FLUX') ;measured flux
mask = mrdfits(modelcube,'MASK') ;mask for flux array
model = mrdfits(modelcube,'MODEL') ;best fit model
emline = mrdfits(modelcube,'EMLINE') ;best fit model with only emlines
emline_base = mrdfits(modelcube,'EMLINE_BASE') ;model of constant baseline from emline fits
;calculate stellar continuum by subtracting emline models from full model
stellar_cont = model - emline - emline_base
;calculate residuals of model
residuals = flux-model
spec_mask = mask[xy[0],xy[1],*]
sel=where(spec_mask eq 0)
cgplot,wave[sel],flux[xy[0],xy[1],sel],xtitle='wavelength (Angstroms)',$
ytitle='flux (1e-17 erg/s/cm^2/Angstrom/spaxel)',/xsty
cgoplot,wave[sel],model[xy[0],xy[1],sel]+5,color='red'
cgoplot,wave[sel],stellar_cont[xy[0],xy[1],sel]+10,color='lime green'
cgoplot,wave[sel],emline[xy[0],xy[1],sel]+15,color='blue'
cgoplot,wave[sel],residuals[xy[0],xy[1],sel]+20,color=cgcolor('orange')
al_legend,['flux','model','stellar continuum','emission lines','residuals'],$
psym=0,color=['black','red','lime green','blue','orange'],/top,/left,charsize=1.5
HYB10 Files
VOR10 files.
file=dir+'manga-7443-12703-MAPS-HYB10-MILESHC-MASTARSSP.fits.gz' snr=mrdfits(file,'BIN_SNR') maxsnr = max(snr,ind) xy = array_indices(snr,ind) ;convert into 2d indices
Recall that although stellar parameters are measured using the Voronoi bins, the emission-line parameters are estimated using the individual spaxels. Therefore, if we want to compare the data to the full best-fitting model that includes stellar continuum and emission lines, we need to use the spectra from the DRP LOGCUBE file, not the DAP HYB10 LOGCUBE file. Let's compare the best-fitting model to the data at the position (j,i) found above (DAP HYB10 LOGCUBE can be downloaded here):
modelcube = dir+'manga-7443-12703-LOGCUBE-HYB10-MILESHC-MASTARSSP.fits.gz' drpcube = dir+'manga-7443-12703-LOGCUBE.fits.gz' flux=mrdfits(drpcube,'flux') fluxmask = mrdfits(drpcube,'MASK') model=mrdfits(modelcube,'model') modelmask=mrdfits(modelcube,'mask') sel=where(fluxmask[xy[0],xy[1],*] eq 0) cgplot,wave[sel],flux[xy[0],xy[1],sel],xtitle='wavelength (Angstroms)',$ ytitle='flux (1e-17 erg/s/cm^2/Angstrom/spaxel)',/xsty sel=where(modelmask[xy[0],xy[1],*] eq 0) cgoplot,wave[sel],model[xy[0],xy[1],sel]+5,color='red' al_legend,['flux','model'],psym=0,color=['black','red'],/top,/left,charsize=1.5
If we just want to compare the best-fitting stellar continuum to the data, we should use the binned spectra within the DAP HYB10 LOGCUBE file:
flux=mrdfits(modelcube,'flux') mask = mrdfits(modelcube,'mask') emline = mrdfits(modelcube,'emline') emline_base=mrdfits(modelcube,'emline_base') stellarcontinuum = model - emline - emline_base sel=where(mask[xy[0],xy[1],*] eq 0) cgplot,wave[sel],flux[xy[0],xy[1],sel],xtitle='wavelength (Angstroms)',$ ytitle='flux (1e-17 erg/s/cm^2/Angstrom/spaxel)',/xsty,yrange=[1.2,2.2] cgoplot,wave[sel],stellarcontinuum[xy[0],xy[1],sel],color='lime green' al_legend,['flux','stellar continuum'],psym=0,color=['black','lime green'],/bottom,/right,charsize=1.5
