Running mangadap on MUSE datacube

[margudo]

Hello,

I'm trying to run the DAP on a MUSE datacube using DataCube class following the provided example for KCWI data but I'm getting a ValueError: No good spectra! when running the SPATIALLY BINNING SPECTRA module.

I'm running dap using the default plan
manga_dap -f ADP-2019-10-08T16-15-37-213.fits --cube_module muse MUSEDataCube --plan_module mangadap.config.analysisplan.AnalysisPlan -o dap

I'm attaching the file with the defined class, in case my error is coming from there.
muse.txt

If someone has trying this already and can point me out to the error I'd appreciate it. Thanks!

[kbwestfall]

Hi @margudo !

That error should indicate that either none of the spaxels met the S/N criteria or they're too highly masked. I looked at your implementation of the MUSEDataCube class, and it looks fine to me. A couple questions for you:

• Did you compare the spaxels in the file to those after the resampling?
• Before the code faulted, did it write the "Reduction QA" reference file? It should be in the common directory with name that's something like *-SNR.fits.gz. You can use that to check the S/N measurements.

I'd like to help get this working. Maybe we can set up a time to chat over zoom, or you can send me the file (via a Google Drive link or something similar) and I can help debug.

[margudo]

Hi Kyle! Spectra looks normal before and after the resampling (besides displacement of emission lines because of the geometric sampling). Looking at the SNR file, 'VARIANCE' is at order 6.40140241e+37 and 'SIGNAL' -9.22337204e+18, so the SNR is negative. I suspect I have a problem with the mask or the units for fluxes. Fluxes are in 10**(-20)*erg/s/cm**2/Angstrom I have tried to not convert to flux densities but it does not solve the problem. I'm masking NaN values but maybe that's too much as you suggest. You can find the original file in the following link: ADP-2019-10-08T16-15-37-213.fits <https://drive.google.com/file/d/14ji_3a5apZGr_rkxp7_nYnPMMb7upEH1/view?usp=drive_web> If it can be helpful we can set up a zoom chat. Thanks, Maria El mié, 6 sept 2023 a las 17:25, Kyle Westfall ***@***.***>) escribió:

…

Hi @margudo <https://github.com/margudo> ! That error should indicate that either none of the spaxels met the S/N criteria or they're too highly masked. I looked at your implementation of the MUSEDataCube class, and it looks fine to me. A couple questions for you: - Did you compare the spaxels in the file to those after the resampling? - Before the code faulted, did it write the "Reduction QA" reference file? It should be in the common directory with name that's something like *-SNR.fits.gz. You can use that to check the S/N measurements. I'd like to help get this working. Maybe we can set up a time to chat over zoom, or you can send me the file (via a Google Drive link or something similar) and I can help debug. — Reply to this email directly, view it on GitHub <#114 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ABZV5PBR3B6V5MLB3UETUOTXZCIYNANCNFSM6AAAAAA4NHI67U> . You are receiving this because you were mentioned.Message ID: ***@***.***>

-- *Maria Argudo-Fernandez, PhD* *Profesora colaboradora* Instituto de Física *|* Pontificia Universidad Católica de Valparaíso *EMERGIA Researcher* Departamento de Física Teórica y del Cosmos Universidad de Granada margudo at ugr.es

[kbwestfall]

Hi Maria,

Thanks for sharing the file! I've been playing around with the (beautiful!) cube this morning. Here's a few things I would try:

• MUSE cubes are very big, so just to test how well the DAP performs, I would only include a few spectra in the output cube. To do this, you can add in the following debug lines immediately after you read the data in the with block:

# DEBUGGING
flux = flux[:,116:119,240:243]
err = err[:,116:119,240:243]
mask = mask[:,116:119,240:243]

• You were defining the mask as a masked array of the flux. This faulted for me, and may be why your VARIANCE and SIGNAL values were so large. Instead, the mask needs to be a boolean array, so I would define mask = numpy.logical_not(numpy.isfinite(flux)).

• The units for MaNGA were 10**(-17)*erg/s/cm**2/Angstrom. Even though I don't think it actually matters in detail, I would convert to these units (divide by 1e3) to put the MUSE data on the same scale (the code probably does assume these units when, e.g., reporting the emission-line fluxes). Also, given that the units are already in flux density (flux per angstrom), you don't need to divide by the wavelength step per pixel (the lines with dw).

• Last, I would reach out to @khrrubin and @erickaguirre96 . They've been working in a branch of the code to enable working with MUSE cubes. I don't know the current status, though.

Here's my edited version of your MUSEDataCube class:

from pathlib import Path
from IPython import embed

import numpy

from astropy.io import fits
from astropy.wcs import WCS

from mangadap.datacube.datacube import DataCube
from mangadap.util.sampling import Resample, angstroms_per_pixel

class MUSEDataCube(DataCube):
    r"""
    Container class for MUSE datacubes.

    Args:
        ifile (:obj:`str`):
            The name of the file to read.
    """

    instrument = 'muse'
    """
    Set the name of the instrument.  This is used to construct the
    output file names.
    """

    def __init__(self, ifile):

        _ifile = Path(ifile).resolve()
        if not _ifile.exists():
            raise FileNotFoundError(f'File does not exist: {_ifile}')

        # Set the paths
        self.directory_path = _ifile.parent
        self.file_name = _ifile.name

        # Collect the metadata into a dictionary
        meta = {}
        meta['z'] = 0.034       # Specific to the target observed
        sres = 2988.98395861011             # Specific to MUSE instrument setting (from header)

        # Open the file and initialize the DataCube base class
        with fits.open(str(_ifile)) as hdu:
            print('Reading MUSE datacube data ...', end='\r')
            prihdr = hdu[1].header
            wcs = WCS(header=prihdr, fix=True)
            flux = hdu[1].data/1e3  # Convert units from 10^-20 (blah) to 10^-17
            err = hdu[2].data/1e3
#            mask = numpy.ma.masked_invalid(flux) # masked nan values
            mask = numpy.logical_not(numpy.isfinite(flux))
        print('Reading MUSE datacube data ... DONE')

        # DEBUGGING
        flux = flux[:,116:119,240:243]
        err = err[:,116:119,240:243]
        mask = mask[:,116:119,240:243]

        # Resample to a geometric sampling
        # - Get the wavelength vector
        spatial_shape = flux.shape[1:][::-1]
        nwave = flux.shape[0]
        coo = numpy.array([numpy.ones(nwave), numpy.ones(nwave), numpy.arange(nwave)+1]).T
        wave = wcs.all_pix2world(coo, 1)[:,2]*wcs.wcs.cunit[2].to('angstrom')

        # - Convert the fluxes to flux density
#        dw = angstroms_per_pixel(wave, regular=False)
#        flux /= dw[:,None,None]
        # - Set the geometric step to the mean value.  This means some
        # pixels will be oversampled and others will be averaged.
        dlogl = numpy.mean(numpy.diff(numpy.log10(wave)))
        # - Resample all the spectra.  Note that the Resample arguments
        # expect the input spectra to be provided in 2D arrays with the
        # last axis as the dispersion axis.
        print('Resampling')
        r = Resample(flux.T.reshape(-1,nwave), e=err.T.reshape(-1,nwave),
                     mask=mask.T.reshape(-1,nwave), x=wave, inLog=False, newRange=wave[[0,-1]],
                     newLog=True, newdx=dlogl)

        # - Reshape and reformat the resampled data in prep for
        # instantiation
        ivar = r.oute.reshape(*spatial_shape,-1)
        mask = r.outf.reshape(*spatial_shape,-1) < 0.8
        ivar[mask] = 0.0
        gpm = numpy.logical_not(mask)
        ivar[gpm] = 1/ivar[gpm]**2
        _sres = numpy.full(ivar.shape, sres, dtype=float)
        flux = r.outy.reshape(*spatial_shape,-1)

        # Default name assumes file names like, e.g., '*_icubew.fits'
        super().__init__(flux, ivar=ivar, mask=mask, sres=_sres,
                         wave=r.outx, meta=meta, prihdr=prihdr, wcs=wcs,
                         name=_ifile.name.split('_')[0])