Getting started

This is your “Quick-Start” guide to working with the DAP output and serves as an introduction to the other, more detailed pages discussing the DAP output. Before diving into the MaNGA data, you are strongly encouraged to briefly familiarize yourself with these pages, in particular.

The DAP provides:

  • Spaxel-by-spaxel coordinate information based on the isophotal ellipticities from the NASA-Sloan Atlas elliptical Petrosian analysis

  • S/N measurements

  • Binned spectra

  • Stellar kinematics and stellar-continuum models

  • Emission-line properties and models

  • Spectral-index measurements

It is important that you understand the procedures and algorithms used by the DAP to provide the data you’re interested in if you are to trust their scientific usage. In particular, consult the DAP Overview paper, Westfall et al. (2019, AJ, 158, 231), the Emission-Line Modeling paper, Belfiore et al. (2019, AJ, 158, 160), and the Data Analysis Pipeline Workflow.

Note

This page is meant to provide a useful guide for how to get started with the DAP output data. If anything is unclear to you then this page is not doing its job!

Please Submit an issue if you have any questions.

Also, if you have helpful code snippets that you’d like to share, please do! The best way to do so would be via a pull request. Please include your code as a script in the examples directory, and include any necessary modules in the mangadap/contrib directory. If necessary, please include licenses relevant to your code in your PR.

Warning

This page is specific to the MaNGA data releases. However, if you’ve used the DAP to analyze a non-MaNGA cube, all of this information is relevant except for the directory structure and the types of analyses performed. Please see How to Fit a non-MaNGA Datacube and Submit an issue if you have problems.


Directory structure

The DAP Directory Structure is described as part of the full DAP Data Model. In particular, note the root directories to the:


DAPTYPE selection

The DAP analyzes each DRP-produced datacube multiple ways. Each DAP Analysis Approach is given a unique keyword, called the DAPTYPE.

The DAPTYPEs in DR15 are:

  • VOR10-GAU-MILESHC: Analysis of spectra binned to \({\rm S/N}\sim10\) using the Voronoi binning algorithm Cappellari & Copin (2003, MNRAS, 342, 345)

  • HYB10-GAU-MILESHC: Stellar-continuum analysis of spectra binned to \({\rm S/N}\sim10\) for the stellar kinematics (same as the VOR10 approach); however, the emission-line measurements are performed on the individual spaxels.

The DAPTYPEs in DR17 are:

  • SPX-MILESHC-MASTARSSP: Analysis of each individual spaxel; spaxels must have a valid continuum fit for an emission-line model to be fit

  • VOR10-MILESHC-MASTARSSP: Analysis of spectra binned to \({\rm S/N}\sim10\) using the Voronoi binning algorithm Cappellari & Copin (2003, MNRAS, 342, 345)

  • HYB10-MILESHC-MASTARSSP: Stellar-continuum analysis of spectra binned to \({\rm S/N}\sim10\) for the stellar kinematics (same as the VOR10 approach); however, the emission-line and spectral-index measurements are performed on the individual spaxels.

  • HYB10-MILESHC-MASTARHC2: Same as the above except the hierarchically clustered MaStar stellar spectra are used to fit the stellar continuum in the emission-line module.

The type of analysis that you should use will depend on your science application. In all cases, please consult Westfall et al. (2019, AJ, 158, 231) and Belfiore et al. (2019, AJ, 158, 160) for usage guidelines and limitations of the data.

SPX-MILESHC-MASTARSSP (DR17 only)

These are useful for most science applications that can push to very low S/N. They are also useful for characterizing the performance of the measurements toward the low S/N limit of the data.

Important considerations:

  • Spectra with \(g\)-band \({\rm S/N} < 1\) will not have a stellar-continuum model or Gaussian emission-line models.

  • Beware of biases in the stellar velocity dispersion with \({\rm S/N}_g < 10\); cf. Westfall et al. (2019, AJ, 158, 231)

VOR10-GAU-MILESHC (DR15), VOR10-MILESHC-MASTARSSP (DR17)

These data are geared toward stellar kinematics, where the data are Voronoi binned to a \(g\)-band \({\rm S/N}\sim 10\).

Important considerations:

  • Spectra with \(g\)-band \({\rm S/N} < 1\) are not included in any bin.

  • Voronoi-binned spectra are just simple means of all the spectra in the bin.

  • The covariance in the datacube is propagated to the variance in the stacked spectra.

  • The spectral resolution in each binned spectra is propagated to the expected spectral resolution in the stacked spectrum, similar to the formalism used by instrumental_dispersion_plane().

  • (Binned) Spectra with \(g\)-band \({\rm S/N} < 1\) will not have a stellar-continuum model or Gaussian emission-line model.

  • Because the binning is done based on the continuum S/N, this limits the emission-line science that can be done at low continuum S/N.

HYB10-GAU-MILESHC (DR15), HYB10-MILESHC-MASTARSSP (DR17), HYB10-MILESHC-MASTARHC2 (DR17)

These are the default files that most users will want to use. We first Voronoi-binned the spectra to a \(g\)-band \({\rm S/N}\sim 10\) to measure the stellar kinematics. Then these bins are deconstructed to fit the emission lines. We refer to this approach as the Hybrid (HYB) binning scheme.

Important considerations:

  • Spectra with \(g\)-band \({\rm S/N} < 1\) are not included in any bin.

  • Voronoi-binned spectra are just simple means of all the spectra in the bin.

  • The covariance in the datacube is propagated to the variance in the stacked spectra.

  • The spectral resolution in each binned spectra is propagated to the expected spectral resolution in the stacked spectrum, similar to the formalism used by instrumental_dispersion_plane().

  • (Binned) Spectra with \(g\)-band \({\rm S/N} < 1\) will not have a stellar-continuum model.

  • All spectra with 80% valid pixels will have a combined emission-line+stellar-continuum model, where the stellar kinematics have been fixed by the fits to the (nearest) binned spectrum.

  • This is the only file where the bin IDs are different for the emission-line properties and spectral indices.


Output files

The primary output files are located at:

SAS Directory

DR15 SAS

DR17 SAS

There are two main output files for each observation (plate-ifudesign combination):

  • manga-[PLATE]-[IFUDESIGN]-MAPS-[DAPTYPE].fits.gz, see DAP MAPS file: 2D “maps” (i.e., images) of DAP measured properties.

  • manga-[PLATE]-[IFUDESIGN]-LOGCUBE-[DAPTYPE].fits.gz, see DAP Model LOGCUBE file: 3D data cubes with the binned and best-fitting-model spectra.

The datacubes produced by the DAP have the same shape as the DRP datacube, and the DAP maps have the same spatial dimensions as a single wavelength channel in the DRP datacubes. This is meant to ease associating the DRP input and DAP output products.

Examples are given below for how to interact with the two main output files using python. However, you are strongly encouraged to install Marvin and use it to interact with the data.


Additionally, the DAP produces a summary file, the DAPall database, with global metrics for each observation and analysis approach (DAPTYPE). This file is at the top-level directory of each MPL/DR.


MAPS files

DAP MAPS file: The MAPS files are the primary output file from the DAP.

In brief, the file contains 2D “maps” (i.e., images) of DAP measured properties. Most properties are provided in groups of three fits extensions:

  1. [property]: the measurement value,

  2. [property]_IVAR: the measurement uncertainty stored as the inverse variance, and

  3. [property]_MASK: a corresponding bit mask for each spaxel.

The headers of each extension provides the astrometric World Coordinate System (WCS) and should exactly match that of the DRP output LOGCUBE files (apart from the wavelength coordinate).

Many properties have multiple “species” or channels associated with them. The identifying name of each mapped property is provided in the header; e.g., the emission lines fit in DR17 are listed here as they appear in the header of the relevant MAPS file extension (e.g., EMLINE_GFLUX). In python, you can create a dictionary of items in each channel using the following method:

# Declare a function that creates a dictionary for the columns in the
# multi-channel extensions
def channel_dictionary(hdu, ext, prefix='C'):
    """
    Construct a dictionary of the channels in a MAPS file.
    """
    channel_dict = {}
    for k, v in hdu[ext].header.items():
        if k[:len(prefix)] == prefix:
            try:
                i = int(k[len(prefix):])-1
            except ValueError:
                continue
            channel_dict[v] = i
    return channel_dict

def channel_units(hdu, ext, prefix='U'):
    """
    Construct an array with the channel units.
    """
    cu = {}
    for k, v in hdu[ext].header.items():
        if k[:len(prefix)] == prefix:
            try:
                i = int(k[len(prefix):])-1
            except ValueError:
                continue
            cu[i] = v.strip()
    channel_units = numpy.empty(max(cu.keys())+1, dtype=object)
    for k, v in cu.items():
        channel_units[k] = v
    return channel_units.astype(str)

This is identical to channel_dictionary(), such that you can do the following:

from mangadap.util.fileio import channel_dictionary
from astropy.io import fits

hdu = fits.open('manga-8138-12704-MAPS-HYB10-MILESHC-MASTARSSP.fits.gz')
emlc = channel_dictionary(hdu, 'EMLINE_GFLUX')
ha_flux = hdu['EMLINE_GFLUX'].data[emlc['Ha-6564']]

In general, it is best to select the extension and channel based on its name, not its number. This is because the ordering of, e.g., the emission lines in the relevant extensions has changed between different DRs/MPLs and may change again. An example of how to do this is shown below.

Usage example

With the MAPS fits file, you should be able to extract DAP maps output using any fits reader. Please Submit an issue if you run into any problems!

For example, we provide a python code snippet below that will plot the \({\rm H}\alpha\) flux map, stellar velocity field, the corrected stellar velocity dispersion field, and the corrected \({\rm H}\beta\) index map for manga-8138-12704-MAPS-HYB10-MILESHC-MASTARSSP.fits.gz. Importantly, note the use of the MAPS Corrections!!

# Imports
import numpy
from matplotlib import pyplot
from astropy.io import fits
from mangadap.util.fileio import channel_dictionary, channel_units

def apply_index_dispersion_correction(indx, indxcorr, unit):
    """
    Apply a set of dispersion corrections.
    """
    if unit not in [ 'ang', 'mag', '' ]:
        raise ValueError('Unit must be mag or ang.')
    return indx * indxcorr if unit in ['ang',''] else indx + indxcorr

# Open the fits file
hdu = fits.open('manga-8138-12704-MAPS-HYB10-MILESHC-MASTARSSP.fits.gz')

# Build dictionaries with the emission-line and spectral-index
# channel names to ease selection and get the spectral-index units
emlc = channel_dictionary(hdu, 'EMLINE_GFLUX')
spic = channel_dictionary(hdu, 'SPECINDEX')
spiu = channel_units(hdu, 'SPECINDEX')

# Show the Gaussian-fitted H-alpha flux map
mask_ext = hdu['EMLINE_GFLUX'].header['QUALDATA']
halpha_flux = numpy.ma.MaskedArray(hdu['EMLINE_GFLUX'].data[emlc['Ha-6564']],
                                   mask=hdu[mask_ext].data[emlc['Ha-6564']] > 0)

pyplot.imshow(halpha_flux, origin='lower', interpolation='nearest', cmap='inferno')
pyplot.colorbar()
pyplot.show()

# Show the stellar velocity field
mask_ext = hdu['STELLAR_VEL'].header['QUALDATA']
stellar_vfield = numpy.ma.MaskedArray(hdu['STELLAR_VEL'].data, mask=hdu[mask_ext].data > 0)

pyplot.imshow(stellar_vfield, origin='lower', interpolation='nearest', vmin=-300, vmax=300,
              cmap='RdBu_r')
pyplot.colorbar()
pyplot.show()

# Show the corrected stellar velocity dispersion field
mask_ext = hdu['STELLAR_SIGMA'].header['QUALDATA']
stellar_sfield_sqr = numpy.ma.MaskedArray(numpy.square(hdu['STELLAR_SIGMA'].data)
                                          - numpy.square(hdu['STELLAR_SIGMACORR'].data[0]),
                                          mask=hdu[mask_ext].data > 0)
# WARNING: This will ignore any data where the correction is larger than the measurement
stellar_sfield = numpy.ma.sqrt(stellar_sfield_sqr)

pyplot.imshow(stellar_sfield, origin='lower', interpolation='nearest', cmap='viridis')
pyplot.colorbar()
pyplot.show()

# Show the corrected H-beta index map
mask_ext = hdu['SPECINDEX'].header['QUALDATA']
hbeta_index_raw = numpy.ma.MaskedArray(hdu['SPECINDEX'].data[spic['Hb']],
                                       mask=hdu[mask_ext].data[spic['Hb']] > 0)
hbeta_index = apply_index_dispersion_correction(hbeta_index_raw,
                                                hdu['SPECINDEX_CORR'].data[spic['Hb']],
                                                spiu[spic['Hb']])

pyplot.imshow(hbeta_index, origin='lower', interpolation='nearest', cmap='inferno')
pyplot.colorbar()
pyplot.show()

Model LOGCUBE files

DAP Model LOGCUBE file: The LOGCUBE files provide the binned spectra and the best-fitting model spectrum for each spectrum that was successfully fit.

These files are useful for detailed assessments of the model parameters because they allow you to return to the spectra and compare the model against the data. As described by Westfall et al. (2019, AJ, 158, 231), the DAP fits the spectra in two stages, one to get the stellar kinematics and the second to determine the emission-line properties. The emission-line module (used for all binning schemes) fits both the stellar continuum and the emission lines at the same time, where the stellar kinematics are fixed by the first fit. The stellar-continuum models from the first fit are provided in the STELLAR extension; to get the stellar continuum determined during the emission-line modeling, you have to subtract the emission-line model (in the EMLINE extension) from the full model (in the MODEL extension); see our Usage example.

Warning

In the HYB binning case the binned spectra provided in the LOGCUBE files are from the Voronoi binning step. However, the emission-line models are fit to the individual spaxels. So:

  • The stellar-continuum fits from the first iteration, in the STELLAR extension, should be compared to the Voronoi binned spectra in the file, but

  • the best-fitting model spectra in the MODEL extension should be compared to the individual spectra from the DRP LOGCUBE file!

Usage example

With the LOGCUBE fits file, you should be able to extract the binned spectra and best-fitting models produced by the DAP using any fits reading software. Please Submit an issue if you run into any problems!

For example, we provide a python code snippet below that plots the highest S/N spectrum, the full model, the residuals, the model stellar continuum, and the model emission-line spectrum using manga-8138-12704-LOGCUBE-HYB10-MILESHC-MASTARSSP.fits.gz. Importantly, note the use of the Maskbits!!

# Imports
import numpy
from astropy.io import fits
from matplotlib import pyplot

# This is a bitmask handling object from the DAP source code
from mangadap.dapfits import DAPCubeBitMask

# Open the fits file
hdu_maps = fits.open('manga-8138-12704-MAPS-SPX-MILESHC-MASTARSSP.fits.gz')
hdu_cube = fits.open('manga-8138-12704-LOGCUBE-SPX-MILESHC-MASTARSSP.fits.gz')

# Get the S/N per bin from the MAPS file
snr = numpy.ma.MaskedArray(hdu_maps['BIN_SNR'].data, mask=hdu_maps['BINID'].data[0] < 0)

# Select the bin/spaxel with the highest S/N and get the relevant map coordiantes
k = numpy.ma.argmax(snr.ravel())
j, i = numpy.unravel_index(k, hdu_maps['BIN_SNR'].data.shape)

# Declare the bitmask object to mask selected pixels
bm = DAPCubeBitMask()
wave = hdu_cube['WAVE'].data
flux = numpy.ma.MaskedArray(hdu_cube['FLUX'].data[:,j,i],
                            mask=bm.flagged(hdu_cube['MASK'].data[:,j,i],
                                                            ['IGNORED', 'FLUXINVALID', 'IVARINVALID',
                                             'ARTIFACT']))

model = numpy.ma.MaskedArray(hdu_cube['MODEL'].data[:,j,i],
                             mask=bm.flagged(hdu_cube['MODEL_MASK'].data[:,j,i], 'FITIGNORED'))
stellarcontinuum = numpy.ma.MaskedArray(hdu_cube['MODEL'].data[:,j,i]
                                            - hdu_cube['EMLINE'].data[:,j,i],
                                        mask=bm.flagged(hdu_cube['MODEL_MASK'].data[:,j,i],
                                                        'FITIGNORED'))
emlines = numpy.ma.MaskedArray(hdu_cube['EMLINE'].data[:,j,i],
                               mask=bm.flagged(hdu_cube['MODEL_MASK'].data[:,j,i], 'ELIGNORED'))
resid = flux-model-0.5

pyplot.step(wave, flux, where='mid', color='k', lw=0.5)
pyplot.plot(wave, model, color='r', lw=1)
pyplot.plot(wave, stellarcontinuum, color='g', lw=1)
pyplot.plot(wave, emlines, color='b', lw=1)
pyplot.step(wave, resid, where='mid', color='0.5', lw=0.5)
pyplot.show()

DAPall file

DAPall database: The DAP provides a single summary file for all analyses done for a given MPL/DR.

The name of the file is dapall-$MANGADRP_VER-$MANGADAP_VER.fits; for example, this is dapall-v3_1_1-3.1.0.fits in DR17.

This file is used primarily for sample selection, and it has many limitations (e.g., it provides star-formation rates calculated without applying an internal extinction correction).

The code below provides an example use, plotting the stellar mass from the DRPall file against the star-formation rates from the DAPall file.

from matplotlib import pyplot
from astropy.io import fits

drpall = fits.open('drpall-v3_1_1.fits')['MANGA'].data
daptype = 'HYB10-MILESHC-MASTARSSP'
dapall = fits.open('dapall-v3_1_1-3.1.0.fits')[daptype]
mass = drpall['nsa_sersic_mass'][dapall['DRPALLINDX']]

indx = (dapall['DAPDONE'] == 1) & (dapall['DAPQUAL'] == 0) & (mass > 0) \
            & (dapall['SFR_TOT'] > -999)

pyplot.scatter(mass[indx], dapall['SFR_TOT'][indx], marker='.', s=20, lw=0, color='k')
pyplot.xlim([1e8, 1e12])
pyplot.ylim([1e-4, 1e1])
pyplot.xscale('log')
pyplot.yscale('log')
pyplot.xlabel(r'$\mathcal{M}_{\ast}$ [$h^{-2} \mathcal{M}_{\odot}$]')
pyplot.ylabel(r'SFR [$h^{-2} \mathcal{M}_{\odot}$/yr]')
pyplot.show()

Using the pixel/spaxel masks

The maskbits for the DAP data are described here. In particular, be aware of the DONOTUSE and the UNRELIABLE flags for the MAPS files. The 2d MAPS file pixel mask is MANGA_DAPPIXMASK. The 3d LOGCUBE file spaxel mask is MANGA_DAPSPECMASK.

In all cases, the DAP has a convenience class that allows a user to quickly determine if any mask value is flagged with a certain value. For example:

# Imports
import os
from astropy.io import fits
from mangadap.util.bitmask import BitMask
from mangadap.config.defaults import sdss_maskbits_file

# Define the path to the IDLUTILS maskbits file
sdssMaskbits = sdss_maskbits_file()

# Instantiate the BitMask object
bm = BitMask.from_par_file(sdssMaskbits, 'MANGA_DAPQUAL')

# Read a DAP file
hdu = fits.open('manga-8138-12704-MAPS-SPX-MILESHC-MASTARSSP.fits.gz')

# Check if the file is critical and print the result
dap_file_is_critical = bm.flagged(hdu['PRIMARY'].header['DAPQUAL'], flag='CRITICAL')
print('This DAP file {0} flagged as CRITICAL.'.format('is' if dap_file_is_critical
                                                      else 'is not'))

There are also a predefined set of derived BitMask classes that the DAP provides. For example:

#Imports
import numpy
from astropy.io import fits
from mangadap.util.drpfits import DRPFitsBitMask

# Instantiate the BitMask object
bm = DRPFitsBitMask()

# Read a DRP file
hdu = fits.open('manga-8138-12704-LOGCUBE.fits.gz')

# Find the number of pixels flagged as DONOTUSE or FORESTAR
indx = bm.flagged(hdu['MASK'].data, flag=['DONOTUSE', 'FORESTAR'])
print('This DRP file has {0}/{1} pixels flagged as either DONOTUSE or FORESTAR.'.format(
        numpy.sum(indx), numpy.prod(indx.shape)))

See also the Marvin Maskbits utility.


Using the BINID extension

The BINID extension has 5 channels. They provide the IDs of spaxels associated with:

  1. each binned spectrum. Any spaxel with BINID=-1 as not included in any bin.

  2. any binned spectrum with an attempted stellar kinematics fit.

  3. any binned spectrum with emission-line moment measurements.

  4. any binned spectrum with an attempted emission-line fit.

  5. any binned spectrum with spectral-index measurements.

In any of these channels, you can obtain the unique bin numbers using numpy.unique(bin_indx.ravel())[1:]; the selection of all but the first array element is just providing all the numbers without the -1 for invalid spaxels (assuming all bin ID maps will have spaxels that are not within the IFU field-of-view, which is always true for MaNGA). If you’re working with anything but the SPX binning, you’ll want to extract the unique spectra and/or maps values. You can do that by finding the indices of the unique bins, like this:

unique_bins, unique_indices = tuple(map(lambda x : x[1:], numpy.unique(bin_indx.ravel(),
                                                                       return_index=True)))

Here’s a worked example where I use unique_bins() to pull out the unique stellar velocities and produce a scatter plot of the x and y positions of the luminosity-weighted bin centers and color them by the measure stellar velocity.

#Imports
import numpy
from astropy.io import fits
from matplotlib import pyplot

from mangadap.util.fitsutil import DAPFitsUtil

# Read a DAP MAPS file
hdu = fits.open('manga-8138-12704-MAPS-HYB10-MILESHC-MASTARSSP.fits.gz')

# Get the unique indices of the stellar-kinematics bins
ubins, uindx = DAPFitsUtil.unique_bins(hdu['BINID'].data[1], return_index=True)

# Get the x and y coordinates and the stellar velocities
x = hdu['BIN_LWSKYCOO'].data[0].ravel()[uindx]
y = hdu['BIN_LWSKYCOO'].data[1].ravel()[uindx]
v = numpy.ma.MaskedArray(hdu['STELLAR_VEL'].data.ravel()[uindx],
                         mask=hdu['STELLAR_VEL_MASK'].data.ravel()[uindx] > 0)

fig = pyplot.figure(figsize=pyplot.figaspect(1))

ax = fig.add_axes([0.15, 0.15, 0.65, 0.65], facecolor='0.8')
cax = fig.add_axes([0.81, 0.15, 0.02, 0.65])
ax.minorticks_on()
ax.set_xlim([18,-18])
ax.set_ylim([-18,18])
ax.grid(True, which='major', color='0.7', zorder=0, linestyle='-')

sp = ax.scatter(x, y, c=v, vmin=-300, vmax=300, cmap='RdBu_r', marker='.', s=30, lw=0, zorder=3)
pyplot.colorbar(sp, cax=cax)

ax.text(0.5, -0.1, r'$\xi$ (arcsec)',
        horizontalalignment='center', verticalalignment='center', transform=ax.transAxes)
ax.text(-0.13, 0.5, r'$\eta$ (arcsec)', rotation='vertical',
        horizontalalignment='center', verticalalignment='center', transform=ax.transAxes)
cax.text(5, 0.5, r'$V_\ast$ (km/s)', rotation='vertical',
         horizontalalignment='center', verticalalignment='center', transform=cax.transAxes)

pyplot.show()