MRK 6: red-shifted HI absorption 5.5.0: Difference between revisions
| Line 223: | Line 223: | ||
default('split') | default('split') | ||
splitms = 'ab658-MKN6.split.ms' | splitms = 'ab658-MKN6.split.ms' | ||
if os.path.exists(splitms): rmtables(splitms) | |||
vis = measurementSet | vis = measurementSet | ||
outputvis = splitms | outputvis = splitms | ||
Revision as of 19:20, 24 February 2010
Overview
This example demonstrates a trickier spectral line data set. In this case, several radio bright active galaxies were observed to look for redshifted 21 cm absorption. The data were obtained using 4IF mode, which means that the full range of velocities were split into two spectral windows (in this example, CASA identifies them as spectral windows 14 and 15). Producing the data cube will require eventually stitching these windows back together. Happily, clean includes flexible regridding and interpolation in the spectral domain. Clean further allows you to adjust the velocity reference frame.
Retrieve the data from the NRAO archive
Download your data from the VLA Archive; in this example we'll use the publicly available survey AB658. These data were published in Gallimore et al. (1999).
With the present archive defaults, you should have obtained the following VLA archive files.
'AB658_A921122.xp1' 'AB658_A921122.xp2' 'AB658_A921122.xp3' 'AB658_A921122.xp4'
This tutorial broadly follows the techniques described in the continuum survey tutorial and the 21 cm emission tutorial, and the basic calibration steps are presented only in outline.
Loading the Data
Recall that the python file globber glob is your friend here!
measurementSet = "ab658.ms" # Store the desired name of the measurement set
from glob import glob
fileList = sorted(glob("AB658*.xp?"))
# double check to see if the file already exists
import os
if not os.path.isdir("./" + measurementSet): # measurement sets are directories
importvla(archivefiles=fileList, vis=measurementSet) # only import if ab658.ms doesn't yet exist
vis = measurementSet
mode = "summary"
vishead
For the purposes of this tutorial, we are interested only in the source Markarian 6, but the reduction techniques could be applied to any of the sources in the measurement set.
Using listobs and plotants we learn the source names, spectral windows, and viable reference antennas. A screenshot of the listobs output is given at right.
| Source | MKN6 |
| Amp Cal | 0134+329 = 3C48 |
| Phase Cal | 1003+830 |
| Central Antennas | VA27, VA09, VA23 |
| Spectral Window IDs | 14, 15 |
It's just as well to put that information into python global variables.
sourceName = "MKN6"
phaseCal = "1003+830"
ampCal = "0134+329"
refAnt = "VA27" # or perhaps "VA09" or "VA23"
spwStr = '14, 15'
The selection of the spectral windows, stored here in the string variable spwStr, is key to getting the calibration right. Each science target in the multisource measurement set has a tuning respective of the source redshift, and the calibration observations are matched to those spectral windows. Were the spectral windows parameter (usually spw) left to its default value (all available spectral windows), many calibration routines would necessarily fail.
Editing the Data
First inspect the data using plotms.
default(plotms)
vis = measurementSet
field = sourceName
plotms()
Be sure to select Axes → X Axis → Time and Axes → Y Axis → Amp. There is some obvious junk isolated in time; see the screenshot at right. Follow the tutorial to flag discrepant data for the source and calibrators.
Calibration
Calibration of HI absorption spectra follows the basic techniques outlined in the continuum tutorial and HI emission tutorial.
First, set the flux of the amplitude calibrator. Notice that we are being careful to keep track of the spectral windows and not using the default spw.
# set the flux standard
default('setjy')
vis = measurementSet
field = ampCal
spw = spwStr # this is key to getting the calibration right!
modimage = ''
setjy()
Bandpass calibration comes next.
# set some default table names
btable = 'ab658.bcal' # bandpass
gtable = 'ab658.gcal' # gain
ftable = sourceName + '.fluxscale' # flux
if os.path.exists(btable): rmtables(btable)
if os.path.exists(gtable): rmtables(gtable)
if os.path.exists(ftable): rmtables(ftable)
# bandpass calibration
default('bandpass')
vis = measurementSet
spw = spwStr # this is key to getting the calibration right!
caltable = btable
gaintable = ''
gainfield = ''
interp = ''
field = phaseCal
spw = ''
selectdata = False
gaincurve = False
opacity = 0.0
bandtype = 'B'
solint = 'inf'
combine = 'scan'
refant = refAnt
bandpass()
In order to perform the amplitude and phase calibrations, we'll need to inspect the bandpass response curve to select appropriate channels to average. We want to choose those channels where the response is relatively flat. The result of plotcal is shown at right: the bandpass is flat over channels 2:27.
# plot the bandpass calibration
default('plotcal')
spw = spwStr
caltable = btable
field = phaseCal
subplot = 211
yaxis = 'amp'
showgui = True
plotcal()
subplot = 212
yaxis = 'phase'
plotcal()
The following steps generate amplitude and phase solutions for the calibrators, and the fluxscale is bootstrapped from the amplitude calibrator onto the phase calibrator, effectively turning the phase calibrator into a local flux calibrator for the science target.
default('gaincal')
vis = measurementSet
caltable = gtable
gaintable = btable
gainfield = ''
interp = 'nearest'
field = ampCal + ',' + phaseCal
# specify data only from the flat region of the bandpass response
spw = '14:2~27,15:2~27'
selectdata = False
gaincurve = False
opacity = 0.0
gaintype = 'G'
solint = 'inf'
combine = ''
calmode = 'ap'
minsnr = 1.0
refant = refAnt
gaincal()
#=====================================================================
#
# Bootstrap flux scale
#
vis = measurementSet
# set the name for the output rescaled caltable
fluxtable = ftable
caltable = gtable
reference = ampCal
transfer = phaseCal
fluxscale()
The task fluxscale will give the bootstrapped flux density of the phase calibrator. If all goes well, something like the following should appear in the log window.
Beginning fluxscale--(MSSelection version)------- Found reference field(s): 0134+329 Found transfer field(s): 1003+830 Flux density for 1003+830 in SpW=14 is: 0.477297 +/- 0.000606617 (SNR = 786.818, nAnt= 26) Flux density for 1003+830 in SpW=15 is: 0.472778 +/- 0.000521238 (SNR = 907.03, nAnt= 26) Storing result in ab658.fluxscale
We should inspect the calibration solutions to make sure they are reasonable, but we'll be cavalier and tread forward. Apply the calibrations to the science target and split that target into its own measurement set.
default('applycal')
vis = measurementSet
gaintable = [ftable,btable]
gainfield = [phaseCal,'*']
interp = ['linear','nearest']
field = sourceName
# choose the appropriate spectral windows
spw = spwStr
# by default, calibrations for a spectral window will map to the same
spwmap = []
selectdata = False
gaincurve = False
opacity = 0.0
applycal()
#
default('split')
splitms = 'ab658-MKN6.split.ms'
if os.path.exists(splitms): rmtables(splitms)
vis = measurementSet
outputvis = splitms
field = sourceName
# weed out the end channels
spw = '14:2~27,15:2~27'
datacolumn = 'corrected'
split()
Continuum Subtraction
Ideally, we would use plotms to identify the channels where HI absorption is present (see the HI emission tutorial.) The absorption line is however weak and isolated against one of a number of radio continuum sources in the data; as such, the absorption is not apparent in averaged visibilities. Instead, we'll have to produce a dirty cube using clean and inspect the spectrum in the image domain.
Search for the Absorption Line
Uvcontsub currently identifies continuum ranges by channel rather than frequency or velocity. We'll make two cubes, one for each spectral window.
default(clean)
vis = splitms
imagename = 'temp0'
spw = '0:0~31'
mode = 'channel'
nchan = -1
start = 0
width = 1
interpolation = 'nearest'
niter = 2000
cell = ['0.4arcsec', '0.4arcsec']
phasecenter = 'J2000 06h52m12.2 +74d25m37'
clean()
#
imagename='temp1'
spw='1:0~31'
clean()
It's not ordinarily necessary to give the phase center, but this is a good opportunity to use the regridding capabilities of clean. The data were originally taken in B1950 coordinates, and the observations followed the good practice of offsetting the pointing center slightly away from the source coordinates. Here, we specify the J2000 coordinates (looked up on NED) and allow clean to handle the precession and centering that we require.
The output cubes will go to temp0.image and temp1.image.
Now we can inspect the cube in viewer.
infile = 'temp0.image'
displaytype = 'raster'
viewer()
Now, use viewer to generate a spectrum. The source spectrum in spw = 1 is shown in the figure at right.
- If necessary. use the player buttons
to scan through the channels until you see the source. - (Tools)Spectral Profile
- Right-click and drag a box around the source; the spectrum will appear in the Image Profile window.
- On the Image Profile window, change the velocity convention to Channel using the pull-down menu.
- Optional: Save the numerical data to a file. Click the save-as-text button
and write the data to a file; for this example, we'll use temp0.txt.
Repeat for spw = 1.
infile = 'temp1.image'
displaytype = 'raster'
viewer()
Perform the Continuum Subtraction
Looking over the plots and text files, only channels 0~25 in either spectral window contain usable data. The absorption line is located in 1:15 and surrounding channels. Just to play it safe, we'll use the channels outside of 1:10~20 (spectral window 1, channels 10-20) to define the continuum for subtraction.
default(uvcontsub)
vis = 'ab658-MKN6.split.ms'
fitspw = '0:0~25,1:0~9;21~25' # define the continuum windows for each spectral window
spw = '0:0~25,1:0~25' # chop out end channels for the continuum and spectral line split
solint = 'int'
fitorder = 1 # fit linear baselines
fitmode = 'subtract'
splitdata = True
uvcontsub()
Uvcontsub produces two new measurement sets:
- ab658-MKN6.split.ms.cont (continuum based on fit to linefree channels), &
- ab658-MKN6.split.ms.contsub (continuum-subtracted visibilities).
Imaging the Continuum Subtracted Cube
Inspecting the preliminary (continuum-containing) cubes, temp0.image and temp1.image, that we generated above, we can determine good parameters for the velocity axis. Here are the velocity options I came up with.
default(clean)
vis = 'ab658-MKN6.split.ms.contsub'
spw=''
imagename = 'MKN6-HIABS'
mode = 'velocity'
nchan = 48
start = '5380km/s'
width = '20km/s'
interpolation = 'linear'
outframe = 'bary'
veltype = 'optical'
It also appears that the channel rms is about 0.5 mJy/beam; to avoid over-cleaning, we'll set the threshold at 1 mJy/beam.
niter = 2000
threshold = '1.0mJy'
weighting = 'natural'
Set the source coordinates. We'll use defaults for the rest of the parameters; with the continuum subtracted, the data have only simple structure with no confusing sources.
# VLA-A; expect 1.5--2 arcsec beams. 0.4arcsec pixels ensure sampling the beam by better than 3x
cell = ['0.4arcsec', '0.4arcsec']
phasecenter = 'J2000 06h52m12.2 +74d25m37'
clean()
The (viewer) results are illustrated in the figure at upper right.
Imaging the Continuum Dataset
We can similarly image the continuum data with just a few changes of the clean parameters.
tget clean
vis = 'ab658-MKN6.split.ms.cont'
mode = 'mfs'
imagename = 'MKN6-CONT'
clean()
The result is shown at right. The continuum image needs some more help. Phase errors produce the 6-pointed star pattern surrounding the source and can be healed using self-calibration. The striping across the field is actually the combination of sidelobes from neighboring sources; these sidelobes can be suppressed by cleaning flanking fields.
Extract the Absorption Spectrum
Viewer provides a handy first-look at the spectrum, but the extracted spectrum is an average over a spatial box. What we really want is the integrated spectrum, not the average. To do this, we'll use CASA's imstat, imval, and imhead tools to extract cube information to python variables and then let python build up the spectrum.
Extract the Spectrum: load the cube data into python arrays
First, we need to determine where in the cube to extract the spectrum; equivalently, we need to find the location of the minimum value in the cube. Imstat is the tool for this.
cubeStat=imstat("MKN6-HIABS.image")
The variable cubeStat is now a python dictionary. To see what it contains, just enter the variable name:
cubeStat
and the following summary of cubeStat appears.
{'blc': array([0, 0, 0, 0], dtype=int32),
'blcf': '06:52:24.903, +74.24.45.777, I, 1.39536e+09Hz',
'flux': array([-0.11339512]),
'max': array([ 0.00302736]),
'maxpos': array([141, 124, 0, 15], dtype=int32),
'maxposf': '06:52:10.909, +74.25.35.400, I, 1.39399e+09Hz',
'mean': array([ -8.62914575e-07]),
'medabsdevmed': array([ 0.00034726]),
'median': array([ -1.18068419e-06]),
'min': array([-0.01676504]),
'minpos': array([127, 131, 0, 9], dtype=int32),
'minposf': '06:52:12.299, +74.25.38.200, I, 1.39454e+09Hz',
'npts': array([ 3145728.]),
'quartile': array([ 0.00069454]),
'rms': array([ 0.00053314]),
'sigma': array([ 0.00053314]),
'sum': array([-2.71449454]),
'sumsq': array([ 0.89414904]),
'trc': array([255, 255, 0, 47], dtype=int32),
'trcf': '06:51:59.574, +74.26.27.778, I, 1.39108e+09Hz'}
These entries can be accessed individually like an array variable. For example, the following commands copy the x and y positions of the absorption line to scalar variables.
x2extract = cubeStat['minpos'][0]
y2extract = cubeStat['minpos'][1]
Now, extract the spectrum using imval.
extractBox = "%d,%d" % (x2extract, y2extract) # copy the position to a string
cubeSpec = imval("MKN6-HIABS.image", box=extractBox, stokes='I')
cubeSpec is another python dictionary, and the spectrum is stored in cubeSpec['data'].
Comment: The absorption line source is unresolved by this observation, or, to put it another way, is entirely contained within a synthesized beam. The extraction at a single spatial pixel gives the integrated spectrum and therefore changes units from surface brightness (mJy/beam) to flux density (mJy). However, the data have been forced into a rectangular grid, and there's no guarantee that the absorption minimum centers within a pixel. Better to use imfit to determine more accurately the position of the minimum, and perhaps re-clean with a new phasecenter. To avoid complicating the present tutorial, we'll just stick with the minimum in the current clean grid.
Generate the frequency axis
We can derive the frequency axis from the header.
import numpy # make sure vector and array arithmetic options are loaded
cubeHead = imhead("MKN6-HIABS.image", mode="list")
nSpec = cubeStat['trc'][3] + 1 # get the number of frequency channels.
f0 = float(cubeHead['crval4']) # reference freq in Hz
df = float(cubeHead['cdelt4']) # channel width in Hz
i0 = cubeHead['crpix4'] # reference pixel
freqSpec = (numpy.arange(nSpec) - i0)*df + f0
Plot the spectrum using matplotlib
Now we have the spectrum stored in freqSpec, cubeSpec['data']. Use matplotlib.pyplot to plot the spectrum (see figure at right).
import matplotlib.pyplot as plt
plt.clf() # clear the plot (figure)
plt.plot(freqSpec/1.e9, cubeSpec['data'], 'k-')
plt.xlabel("Freq (GHz)")
plt.ylabel("Flux Density (mJy)")
A publication-quality figure: plot the spectrum against velocity
Nice, now let's make a velocity axis instead. We'll use the optical convention to convert frequency to velocity; the radio convention is offered as an alternative but is commented out.
fRest = cubeHead['restfreq'][0]
velocity = 299792.458 * (fRest/freqSpec - 1.0) # optical convention, km/s
# velocity = 299792.458 * (1.0 - freqSpec/fRest) # radio convention, km/s
Plot the spectrum again, but this time try for publication-ready quality (see figure at right). This example is somewhat more clever to highlight the use of TeX typography in matplotlib.
# set up fonts
from matplotlib import rc
rc('font',**{'family':'serif','serif':['Times']})
rc('text', usetex=True)
plt.clf() # clear the plot (figure)
plt.plot(velocity, cubeSpec['data'], 'k-')
# the r prefix to the following format strings prevents \ from triggering an escape sequence,
# like \n for line feed
plt.xlabel(r"$cz$ (km s$^{-1}$)",fontsize=16)
plt.ylabel(r"$S_{\nu}$ (mJy)",fontsize=16)
plt.ylim(-0.023,0.005)
plt.savefig("mkn6-abs.eps") # save the figure as encapsulated postscript, suitable for journals
plt.savefig("mkn6-abs.png") # save the figure as a bitmap, better suited for astro-ph
Save the spectrum to a text file
It's handy to have a simple, columnar text file to inspect and analyze later. Or perhaps you might want to analyze the spectrum outside of CASA. Numpy makes it easy to save and restore text files.
# First, combine the spectrum into a single array
spectrum = numpy.vstack((velocity, cubeSpec['data']))
numpy.savetxt("mrk6.txt", numpy.transpose(spectrum))
Here's a listing of the first few lines of mrk6.txt.
5.379999999489542461e+03 -5.940478877164423466e-04 5.399939886933506386e+03 -3.919676819350570440e-04 5.419882380281690530e+03 -3.093948180321604013e-04 5.439827480044931690e+03 -1.092202146537601948e-03 5.459775186734200361e+03 -4.310846794396638870e-03 5.479725500860732609e+03 -3.458642400801181793e-03 5.499678422935699928e+03 -7.933050510473549366e-04 5.519633953470603956e+03 -2.004600130021572113e-03 5.539592092976883578e+03 -8.163747377693653107e-03 5.559552841966238702e+03 -1.676503755152225494e-02
Here's how to restore that text file back into CASA (or python, for that matter).
# import numpy # if you haven't already
spectrum = numpy.loadtxt("mrk6.txt",dtype="float")
x = spectrum[:,0]
y = spectrum[:,1]
--Jack Gallimore 00:36, 12 February 2010 (UTC)
