Stacking Multiple Spectral Lines at Same Position

↵ Hints, Tips, and Tricks

This guide presents a method in CASA for averaging together multiple transitions of the same line. This process increases the signal to noise of the line observations. It may not be suitable for all science cases.

The observations in this guide were of 5 different hydrogen radio recombination lines. The correlator was configured to observe one spectral line per spectral window. The data were calibrated and split off so that each spectral line had its own visibility set.

The lines are all emitted from the same region in the galaxy. Therefore, a strategy to combine them is to use clean to put each spectral line cube on the same velocity scale. Then the cubes are smoothed to the same resolution as the cube with the coarsest resolution using imsmooth. Finally, the spectral line cubes are averaged together using immath.

This script assumes that you have already continuum subtracted the data using uvcontsub. If you need to continuum subtract in the image domain using imcontsub, it should be done after the clean, not after the spectral line cubes are averaged together. Be careful with your continuum subtraction -- small errors in continuum subtraction can lead to spurious lines.

Instead of smoothing after cleaning, another option would be to have all the images use the same restoring beam size. You can set this parameter in clean.

Below we do a detailed walkthrough of the script.

First, we define which spectral line visibilities we are going to use and their rest frequencies.

linems = ['ic342_cband_line_if2_cal_spw0.ms.contsub',
          'ic342_cband_line_if2_cal_spw1.ms.contsub',
          'ic342_cband_line_if2_cal_spw2.ms.contsub',
          'ic342_cband_line_if2_cal_spw3.ms.contsub',
          'ic342_cband_line_if2_cal_spw4.ms.contsub']
restfreq = ['6.2891GHz',
            '6.4788GHz',
            '6.6761GHz',
            '6.8815GHz',
            '7.0954GHz']

Now we define the parameters that we are going to use to image the spectral line data. Depending on your data set, you may want to change some of these parameters.

velstart = '-250km/s'     # velocity to start the cube at
chanwidth = '10km/s'      # width of channels in the cube
nchan = 50                # number of channels in the cube
outframe = 'bary'         # use a barycentric reference frame
veltype = 'optical'       # use the optical velocity definition
imsize = [1536,1536]
cell = ['0.75arcsec']

Now set the parameters for a shallow clean on the data.

threshold = '20mJy'
niter = 300

Finally, we set the name of the output image

output_stacked_image = 'RRL_if2_stack.image'

Now we run the script.

# creating the output image names
output_line_image = [linems[i].replace('.ms.contsub','.contsub') for i in range(len(linems))]

# creating the images for each spectral line data set
for i in range(len(linems)):
    
    clean(vis=linems[i],
          imagename=output_line_image[i],
          mode='velocity',
          nchan=nchan,
          start=velstart,
          width=chanwidth,
          outframe=outframe,
          veltype=veltype,
          niter=niter,
          threshold=threshold,
          interactive=False,
          imsize=imsize,
          cell=cell,
         restfreq=restfreq[i])

# Getting the sizes of the beam.
beammajor = []
beamminor = []
beampa = []

for image in output_line_image :
    myimagename = image + '.image'
    beammajor.append(imhead(imagename=myimagename,mode='get',hdkey='beammajor'))
    beamminor.append(imhead(imagename=myimagename,mode='get',hdkey='beamminor'))
    beampa.append(imhead(imagename=myimagename,mode='get',hdkey='beampa'))

# determining what size the spectral line images should be smoothed to.
targetbeammajor = max(beammajor)
targetidx = beammajor.index(targetbeammajor)
targetbeamminor = beamminor[targetidx]
targetbeampa = beampa[targetidx]

# smoothing the images
for image in output_line_image:
    inputimage = image + '.image'
    outputimage = image + '.image' + '.smooth'
    if image != output_line_image[targetidx]:
        imsmooth(imagename=inputimage,outfile=outputimage,kernel='gauss',
                 major=str(targetbeammajor['value'])+targetbeammajor['unit'],
                 minor=str(targetbeamminor['value'])+targetbeamminor['unit'],
                 pa=str(targetbeampa['value'])+targetbeampa['unit'],
                 targetres=True)
    else:
        # need to skip the one with the largest resolution or imsmooth fails.
        os.system('cp -ir '+ inputimage + ' ' + outputimage)


# creating the expression for immath. It doesn't appear to know about mean.
output_line_image = [output_line_image[i] + '.image' for i in range(len(output_line_image))]
myim = ['IM' + str(i) for i in range(len(output_line_image))]
myexp =  '(' + '+'.join(myim) + ')/' + str(float(len(output_line_image)))

# combining the images
immath(imagename=output_line_image,outfile=output_stacked_image,mode='evalexpr',expr=myexp)

↵ CASAguides