Antennae Band7 - Imaging for CASA 3.3

Overview

This section of the AntennaeBand7 CASA Guide cover the imaging of the continuum and spectral line data. It begins where the Antennae Band7 - Calibration section left off. If you completed the Calibration section of the guide, then you may want to continue with the calibrated data sets. If you did not complete the Calibration section, then you can download the calibrated uv-data from the region closest to your location:

North America Europe East Asia

Then download the file 'Antennae_Band7_CalibratedData.tgz' to obtain the calibrated uvdata. Once the download has finished, unpack the file using 'tar -xvzf Antennae_Band7_CalibratedData.tgz' and change directory (cd) to Antennae_Band7_CalibratedData. You should have Antennae_South.ms and Antennae_North.ms in your working directory. Then start CASA with the 'casapy' command.

Continuum Imaging

Fig. 1. Southern mosaic: Amplitude vs. channel. The CO(3-2) line is seen from 50 to 100

Fig. 2. Northern mosaic: Amplitude vs. channel. The CO(3-2) line is seen from 70 to 100

We will make 345 GHz continuum images for the two regions covered by the mosaics. We use the task clean over the channels that are free of the line emission; we avoid the edge channels which tend to be noisier due to bandpass rolloff effects. The line-free channels are found by plotting the average spectrum. We find the CO(3-2) line from channels 50 to 100 in the southern mosaic, and from 70 to 100 in the northern mosaic.

# In CASA
plotms(vis='Antennae_South.cal.ms',xaxis='channel',yaxis='amp',
      avgtime='1e8',avgscan=T)
plotms(vis='Antennae_North.cal.ms',xaxis='channel',yaxis='amp',
      avgtime='1e8',avgscan=T)

The *avgtime* is set to a large value so that it averages over all the integrations, and *avgscan* is set to allow averaging of the different scans.

Next we create continuum images from the line-free channels.

Northern Continuum Mosaic

For illustrative purposes we first make a dirty image to see if there is emission and what the exact beam size is. It should be on the order of 1" but this will vary a bit according to the uv-coverage in the actual data. We will start with a cell size of 0.2" to oversample the beam by a factor of 5. The imsize needs to be large enough given the cell size to comfortably encompass the mosaic. From the mosaic footprints shown in the overview we can see that Northern mosaic imsize needs to be about 1 arcmin. At 2" this requires a imsize=300.

For mosaics it is essential to pick the center of the field to be imaged using the phasecenter parameter. Otherwise it will default to the first pointing included in the imaging -- since this is often at one corner of the mosaic, the image will not be centered. For the Norther mosaic, the center pointing corresponds to field id=12. You can also set an explicit coordinate (see the clean help for syntax.

Other essential clean parameters for this case include:

• vis='Antennae_North.cal.ms': The calibrated dataset on the science target
• imagename='Antennae_North.Cont.Clean': The base name of the output images:
  • <imagename>.image # the final restored image
  • <imagename>.flux # the effective primary beam response (e.g. for pbcor)
  • <imagename>.flux.pbcoverage # the primary beam coverage (ftmachine=’mosaic’ only)
  • <imagename>.model # the model image
  • <imagename>.residual # the residual image
  • <imagename>.psf # the synthesized (dirty) beam
• spw='0:1~50;120~164': To specify only the line-free channels of spectral window 0.
• mode='mfs': Multi-Frequency Synthesis: The default mode, which produces one image from all the specified data combined. This will grid each channel independently before imaging. For wide bandwiths this will give a far superior result compared to averaging the channels and then imaging.
• niter=0: Maximum number of clean iterations -- niter=0 will do no cleaning

# In CASA
os.system('rm -rf Antennae_North.Cont.Dirty*')
clean(vis='Antennae_North.cal.ms',imagename='Antennae_North.Cont.Dirty',
     field='',phasecenter='12',
     mode='mfs',restfreq='345.79599GHz',
     spw='0:1~50;120~164',
     imagermode='mosaic',
     imsize=300,cell='0.2arcsec',
     interactive=F,niter=0)

The reported beam size is 1.26" x 0.66" P.A. 86.045 degrees.

# In CASA
viewer('Antennae_North.Cont.Dirty.image')

Fig. 3. Residual for Northern continuum mosaic after 100 iterations; the clean mask is shown by the white contour.

Yes there is definitely a detection in the vicinity of the Northern nucleus (see Fig. 4 below). Using the square region icon, and drawing a box near but not including the emission, we find the rms noise is about 0.4 mJy/beam in the dirty image.

Next we switch to refined values of cell=0.13" and imsize=500 based on the observed beam size. We also switch to interactive mode so that you can create a clean mask using the polygon tool (note you need to double click inside the polygon region to activate the mask. See here for a more complete description of interactive clean and mask creation.

• niter=1000: Maximum number of clean iterations -- we will stop interactively
• threshold='0.4mJy': Stop cleaning if the maximum residual is below this value (the dirty rms noise)
• interactive=T: Clean will be periodically interupted to show the residual clean image. Interactive clean mask can be made. If no mask is created, no cleaning is done.
• npercycle=50: How many iterations to clean before returning to interactive clean window.

# In CASA
os.system('rm -rf Antennae_North.Cont.Clean*')
clean(vis='Antennae_North.cal.ms',imagename='Antennae_North.Cont.Clean',
     field='',phasecenter='12',
     mode='mfs',restfreq='345.79599GHz',
     spw='0:1~50;120~164',
     imagermode='mosaic',
     imsize=500,cell='0.13arcsec',
     interactive=T,npercycle=50,
     niter=1000, threshold='0.4mJy')

The residuals are "noise-like" after only 50 iterations so hit the red X symbol in the interactive window to stop cleaning here.

Note that if you run the clean task again with the same imagename, without deleting the existing <imagename>.* files, clean assumes that you want to continue cleaning the existing images. We put an rm command before the clean command to guard against this, but sometimes this is what you want.

Southern Continuum Mosaic

Fig. 3. Residual for Southern continuum mosaic after 100 iterations; the clean mask is shown by the white contour.

For the southern mosaic we modify the phase center, mosaic size, and the line-free channels to be consistent for this mosaic. As above, the mosaic size can be judged from here. We expect the beam to be about the same and use the same cell size. During clean repeat process of making a clean mask.

# In CASA
os.system('rm -rf Antennae_South.Cont.Clean*')
clean(vis='Antennae_South.cal.ms',imagename='Antennae_South.Cont.Clean',
     field='',phasecenter='15',
     mode='mfs',restfreq='345.79599GHz',
     spw='0:1~30;120~164',
     imagermode='mosaic',
     imsize=750,cell='0.13arcsec',
     interactive=T,npercycle=50,
     niter=1000, threshold='0.4mJy')

The beam size reported in the logger for the Southern mosaic: 1.15"x 0.66", and P.A.=68 deg; a little better than the for the North due to better uv-coverage.

Stop after 50 iterations.

Fig. 5. 345 GHz continuum image of the northern mosaic

Fig. 6. 345 GHz continuum image of the southern mosaic

Continuum Image Analysis

We will determine some statistics for the images using the task imstat:

# In CASA
imstat('Antennae_North.Cont.Clean.image')
imstat('Antennae_South.Cont.Clean.image')

The peak flux densities of the images are ~4.4 and 2.9 mJy for the Southern and Northern mosaics, respectively, and both have a similar rms ~ 0.4 mJy, thus the obtained dynamic ranges are larger or equal than 10. Note that if the emission had been strong we should run twice, once excluding the emission from the rms calculation and once with it included to get the peaks using the box command. It is often easier to use either imview or viewer to display the image and then interactively use the region tools to get the statistics.

Next make hardcopy plots of the continuum images using imview (a task wrapper for the viewer to obtain Figures 4 and 5:

# In CASA
imview (
  raster={'file': 'Antennae_North.Cont.Clean.image',
  'colorwedge':T,},
  zoom=1, out='Antennae_North.Cont.Clean.image.png')
imview (
  raster={'file': 'Antennae_South.Cont.Clean.image',
  'colorwedge':T,},
  zoom=1, out='Antennae_South.Cont.Clean.image.png')

CO(3-2) Imaging

Continuum subtraction

In these data, the continuum emission is too weak to contaminate the line emission (i.e. the peak continuum emission is less than the rms noise in the spectral line channels). For illustrative purposes we demonstrat how to subtract the continuum emission in the uv-domain using the task uvcontsub2.

# In CASA
uvcontsub2(vis = 'Antennae_North.cal.ms', fitspw='0:1~50;120~164', fitorder = 1)

uvcontsub2(vis = 'Antennae_South.cal.ms', fitspw='0:1~30;120~164', fitorder = 1)

where: fitspw gives the line-free channels for each mosaic and fitorder=1. Higher order fits are not recommended. If you do not have line-free channels on both sides of the line fitorder=0 is recommended.

Clean CO(3-2) line

Northern CO(3-2) uv-spectrum in LSRK velocity space.

Southern CO(3-2) uv-spectrum in LSRK velocity space.

Now we are ready to make cubes of the line emission. The imaging parameters are similar to the continuum except for those dealing with the spectral setup: mode, start, width, nchan, restfreq, and outframe parameters. When making spectral images you have three choices for the mode parameter: channel, velocity, and frequency. Data are taken using constant frequency channels. For spectral line analysis it's often more useful to have constant velocity channels, and this is also the best way to make images of multiple lines with the exact same channelization for later comparison. For mode='velocity', the desired start and width also need to be given in velocity units for the desired output frame.

It is important to note that ALMA does not do on-line Doppler Tracking and the native frame of the data is TOPO. If you do not specify outframe the output cube will also be in TOPO, which is not very useful for spectral line work. The Doppler Shift is taken out during the regridding to the desired outframe in clean or alternatively it can be done separately by the cvel task which would need to be run before clean.

To see what velocity parameters you want to set in clean it is useful to make a plot in plotms in the desired output frame. Note these plots take a little longer because of the frame shift. In order to compare with recently SMA data, we chose LSRK, but it should be noted that many papers of this source are in the BARY frame.

# In CASA
plotms(vis='Antennae_North.cal.ms.contsub/',xaxis='velocity',yaxis='amp',
       avgtime='1e8',avgscan=T,transform=T,freqframe='LSRK',
       restfreq='345.79599GHz',plotfile='North_CO3_2_vel.png')

# In CASA
plotms(vis='Antennae_South.cal.ms.contsub',xaxis='velocity',yaxis='amp',
      avgtime='1e8',avgscan=T,transform=T,freqframe='LSRK',
      restfreq='345.79599GHz',plotfile='South_CO3_2_vel.png')

As before, it is very important that you make a clean mask. There are many ways to do this ranging from the complicated to simple. For this example we will make a single clean mask that encompasses the line emission in every channel and apply it to all channels. This is much better than no clean mask, though not quite as good as making an individual mask for each channel.

Notable parameters include:

• mode='velocity',outframe='LSRK',restfreq='345.79599GHz'. We use 'velocity' mode to set velocity information in the Local Standard of Rest frame, and use the rest frequency of the CO(3-2) line.
• nchan=60,start='1300km/s',width='10km/s': To produce a data cube with 60 channels("nchan"=60), starting at 1300km/s and with velocity widths of 10 km/s. That will include all the CO(3-2) emission in both mosaics.
• imsize=1000,cell='0.13arcsec'. An image size of 130 arcsec, with pixels of 0.13 arcsec (~ 5 times the minor axis of the synthesized beam). We make the imsize larger than the mosaic to increase the dirty beam image size, as the emission is quite extended along the mosaic.
• wighting='briggs',robust=0.5: Weighting to apply to visibilities. We use Brigg's weighting and robustness parameter 0.5 (in between natural and uniform weighting).
• niter=10000: Maximum number of clean iterations.
• threshold='5.0mJy': Stop cleaning if the maximum residual is below this value.

The threshold is ~ 2 the rms of the line-free channel for each dataset.

Fig. X. Interactive cleaning of northern mosaic. The clean mask is shown in white. In this example we used the same mask for all channels by selecting "all channels" before drawing mask.

# In CASA
os.system('rm -rf Antennae_North.CO3_2Line*')
clean(vis='Antennae_North.cal.ms',
     imagename='Antennae_North.CO3_2Line.Clean',
     spw='0',field='',phasecenter='12',
     mode='velocity',outframe='LSRK',restfreq='345.79599GHz',
     nchan=70,start='1200km/s',width='10km/s',
     imagermode='mosaic',
     imsize=500,cell='0.13arcsec',minpb=0.2,
     interactive=T,
     weighting='briggs',robust=0.5,
     niter=10000, threshold='5.0mJy')

The figure shows what the final clean box should look like. Now clean using the green circle arrow. You can watch the progress in the logger. When 100 iterations are done, the viewer will show the residual map for each channel. Cycle through the channels and see whether you are still happy with the clean box in each signal channel. The "erase button" can help you fix mistakes. If necessary adjust. It is often useful to adjust the colorscale with the "plus" symbol icon. To make it go a bit faster you can increase iteration interactively, to 200 or 300, but don't overdo it.

These data are pretty severely dynamic range limited. In other words the noise in bright channels is set by a maximum signal-to-noise. Also, we are going to self-calibrate the data so its best to be conservative here - it can be difficult in images like these to discern real features from artifacts. If in doubt, don't include it in the clean mask. Keep cleaning and see if the feature becomes weaker you cannot get rid of real emission by not masking it. You can create a brighter feature by masking an artifact.

Stop cleaning after about 800 iterations (Red x).

Inspect the resulting data cube:

# In CASA
imview (raster={'file': 'Antennae_North.Line.Clean.image'})

Repeat process for Southern mosaic.

# In CASA
os.system('rm -rf Antennae_South.CO3_2Line*')
clean(vis='Antennae_South.cal.ms',
     imagename='Antennae_South.CO3_2Line.Clean',
     spw='0',field='',phasecenter='15',
     mode='velocity',outframe='LSRK',restfreq='345.79599GHz',
     nchan=70,start='1200km/s',width='10km/s',
     imagermode='mosaic',
     imsize=750,cell='0.13arcsec',minpb=0.2,
     interactive=T,
     weighting='briggs',robust=0.5,
     niter=10000, threshold='5.0mJy')

Inspect the resulting data cube:

# In CASA
imview (raster={'file': 'Antennae_South.Line.Clean.image'})

The (line-free) channel rms for the northern and southern mosaics are 2.7 mJy/beam and 2.4 mJy/beam, respectively, for 2.3 and 3.8 hours times on source approximately. The ALMA sensitivity calculator gives an on-source observing time of ~3-4 hours to reach a noise level of 2.5 mJy using 12 antennas for an angular resolution of 1 arcsec. Note that the fields are overlapping by 0.5 the size of the primary beam. Note however that the dynamic range is worse than that. Self-calibration and multiscale will help improving this situation.

Self calibration and multi scale cleaning

We will update this section soon.

Moment maps

Next, we make the integrated intensity maps and velocity fields of the CO(3-2) emission using the task immoments. Note that we choose to run the immoments separately for the total intensity map and the velocity field because we want to use different flux thresholds.

For the northern mosaic the moment 0 and 1 can be calculated as follows:

Fig. 6. The CO(3-2) total intensity map (moment 0) of the northern mosaic

Fig. 7. The CO(3-2) velocity field (moment 1) of the northern mosaic

# In CASA
os.system('rm -rf Antennae_North.Line.Clean.Mom.image*') 
immoments('Antennae_North.Line.Clean.image', 
   moments=[0], axis="spec", 
　outfile='Antennae_North.Line.Clean.Mom.image.integrated')

os.system('rm -rf Antennae_North.Line.Clean.Mom.Mask.image.weighted_coord*') 
immoments('Antennae_North.Line.Clean.image', 
   moments=[1], axis="spec", 
    mask="Antennae_North.Line.Clean.image>4.e-2",
   outfile='Antennae_North.Line.Clean.Mom.Mask.image.weighted_coord')

And for the southern mosaic:

Fig. 8. The CO(3-2) total intensity map (moment 0) of the southern mosaic

Fig. 9. The CO(3-2) velocity field (moment 1) of the southern mosaic

# In CASA
os.system('rm -rf Antennae_South.Line.Clean.Mom.image*') 

immoments('Antennae_South.Line.Clean.image', 
   moments=[0], axis="spec",
   outfile='Antennae_South.Line.Clean.Mom.image.integrated')

immoments('Antennae_South.Line.Clean.image', 
   moments=[1], axis="spec", 
   mask="Antennae_South.Line.Clean.image>2.5e-2",
   outfile='Antennae_South.Line.Clean.Mom.Mask.image.weighted_coord')

where:

• moments=[x]: To specify that we wish to make the xth moment map. The 0th moment map gives integrated emission and the 1st gives the intensity-weighted velocity field
• axis='spectral': Indicates the moment axis; in this case, 'spectral'
• mask="Antennae_North.Line.Clean.image>X" Mask applied to the image before calculating the moments. Consider pixels in Antennae_North.Line.Clean.image that are larger than X
• outfile='Antennae_South.Line.Clean.Mom.image.XX': The output image name

Finally, we visualize the different moments using imview:

# In CASA
imview (
  raster={'file': 'Antennae_North.Line.Clean.Mom.image.integrated',
  'colorwedge':T,},
  zoom=1, out='Antennae_North.Line.Clean.Mom.image.integrated.png')

imview (
   raster={'file': 'Antennae_North.Line.Clean.Mom.Mask.image.weighted_coord',
   'range': [1500,1750],'colorwedge':T},
   zoom=1, out='Antennae_North.Line.Clean.Mom.Mask.image.velfield.png')

imview (
   raster={'file': 'Antennae_South.Line.Clean.Mom.image.integrated',
   'colorwedge':T},
   zoom=1,out='Antennae_South.Line.Clean.Mom.image.integrated.png')

imview (
   raster={'file': 'Antennae_South.Line.Clean.Mom.Mask.image.weighted_coord',
   'range': [1300,1900],'colorwedge':T},
   zoom=1, out='Antennae_South.Line.Clean.Mom.Mask.image.velfield.png')

Comparison with previous SMA CO(3-2) data

We compare with SMA CO(3-2) data (Ueda, Iono, Petitpas et al., in preparation, to be submitted to ApJ). Figure 13 shows a comparison plot betwee the moment 0 maps of ALMA and SMA data using the viewer. The fluxes, peak locations and large scale structure are consistent. Both southern and northern components have been combined.

Exporting to Fits

You can convert from CASA format to FITS if you want to analyze the data using another software package. If you have difficulty with the standard FITS output from CASA, try setting velocity=T and/or dropstokes=T.

# In CASA
exportfits(imagename='Antennae_North.Line.Clean.image',
     fitsimage='Antennae_North.Line.Clean.image.fits')
exportfits(imagename='Antennae_North.Line.Clean.Mom.image.integrated',
     fitsimage='Antennae_North.Line.Clean.Mom.image.integrated.fits')
exportfits(imagename='Antennae_North.Line.Clean.Mom.Mask.image.weighted_coord',
     fitsimage='Antennae_North.Line.Clean.Mom.Mask.image.weighted_coord.fits')

# In CASA
exportfits(imagename='Antennae_South.Line.Clean.image',
     fitsimage='Antennae_South.Line.Clean.image.fits')
exportfits(imagename='Antennae_South.Line.Clean.Mom.image.integrated',
     fitsimage='Antennae_South.Line.Clean.Mom.image.integrated.fits')
exportfits(imagename='Antennae_South.Line.Clean.Mom.Mask.image.weighted_coord',
     fitsimage='Antennae_South.Line.Clean.Mom.Mask.image.weighted_coord.fits')