`max_map` calculated from `cutoutcube` but used with `noisemap` that is calculated using `cube`

[jmangum]

max_map is derived from a masked version of cutoutcube, which is the so-called "brightest line cube":

    cutoutVcube = cutoutcube.with_spectral_unit(u.km/u.s,
                                                   rest_value=brightest_line_frequency,
                                                   velocity_convention='optical')

    # Use the brightest line to identify the appropriate peak velocities, but ONLY
    # from a slab including +/- width:
    brightest_cube = cutoutVcube.spectral_slab(vz-velocity_half_range,
                                               vz+velocity_half_range)

    # compute various moments & statistics along the spcetral dimension
    peak_velocity = brightest_cube.spectral_axis[brightest_cube.argmax(axis=0)]
    max_map = peak_amplitude = brightest_cube.max(axis=0) # This is sometimes as all-NaN slice

max_map is used later in cubelinemoment_multiline to calculate peak_sn = max_map / noisemap. noisemap is derived from cube(the target line cube), not from cutoutcube, in cubelinemoment_setup:

noisemap = cube.with_mask(mask[:,None,None]).std(axis=0)

peak_sn = max_map / noisemap

This looks inconsistent...

[keflavich]

Right, this can be inconsistent if the cubes are produced with different channel widths or from a different data set.

[jmangum]

I think that the current behaviour was ok when cube and cutoutcube were always the same. Now that we allow for cube and cutoutcube to be different we should use max_map calculated from cube when calculating peak_sn. Note, though, that there are two options. Since peak_sn = max_map / noisemap, either max_map needs to be changed to use cube or noisemap needs to be changed to use cutoutcube. Since peak_sn is used for setting a threshold for the gaussian mask, I think that changing max_map to use cube appears to be the right choice.

[jmangum]

I have been doing quite a few tests this morning trying to make the peak_sn = max_map / noisemap consistent by using information from just cube or cutoutcube to calculate peak_sn. In the end found that basing calculation on cube resulted in the width_mask, which has a threshold calculated as 1 / peak_sn, was far to restrictive, basically excluding all emission toward most sample pixels. I think that I was not thinking about how peak_sn was being used. Suggest that using max_map from cutoutcube and noisemap from cube is correct for setting the threshold for a theoretical gaussian mask for a hypothetical spectral line that is based on the "tracer" line from cutoutcube. If @keflavich agrees, we can close this issue.

[keflavich]

I think that sounds right.

It's sounding more and more like we need to formally write up this process as a software paper or something, and show figures detailing each step. Though, this writeup is similar: https://github.com/abulatek/tutorials/blob/new/SpectralCubeReprojectMaskExample_v4p1.ipynb

[jmangum]

Agreed. Created Issue #42 to track.