M100 Band3 SingleDish 4.2.2
M100 Single Dish Data Reduction (under modification by AH)
Overview
This portion of the M100 Single Dish Data Reduction CASA Guide will cover the reduction of the Total Power (TP) array data into units of Kelvins on the antenna temperature (Ta*) scale and imaging. Converting this image to the Jansky scale (Jy/beam) to be combined with interferometric data is covered in the M100 Band3 Combine 4.1 section.
This guide is designed for CASA 4.1.0.
If you haven't downloaded the data, you can XXXXXXX and XXXXXX:
Once the download has finished, upack the file:
# In a terminal outside CASA
tar -xvzf XXXsingledish_datasetXXX_TBD.tgz
cd XXrelevant_directoryXXX
# Start CASA
casapy
Summary of the observing
There were six observations made. The table below indicates the uid reference of each and the start and end times.
uid___A002_X60b415_X39a Observed from 14-Apr-2013/05:34:15.7 to 14-Apr-2013/05:58:23.8 (UTC) uid___A002_X60b415_X6f7 Observed from 14-Apr-2013/06:23:03.0 to 14-Apr-2013/06:47:11.0 (UTC) uid___A002_X6218fb_X264 Observed from 28-Apr-2013/04:12:06.1 to 28-Apr-2013/04:36:07.5 (UTC) uid___A002_X6218fb_X425 Observed from 28-Apr-2013/04:38:56.8 to 28-Apr-2013/05:03:00.3 (UTC) uid___A002_X6321c5_X3a7 Observed from 12-May-2013/02:22:16.9 to 12-May-2013/02:43:59.8 (UTC) uid___A002_X6321c5_X5ca Observed from 12-May-2013/02:47:16.8 to 12-May-2013/03:09:00.9 (UTC)
Which version of CASA to use
This guide has been written for CASA release 4.1.0. Please confirm your version before proceeding.
# In CASA
version = casadef.casa_version
print "You are using " + version
if (version < '4.1.0'):
print "YOUR VERSION OF CASA IS TOO OLD FOR THIS GUIDE."
print "PLEASE UPDATE IT BEFORE PROCEEDING."
else:
print "Your version of CASA is appropriate for this guide."
Initial Inspection
We will eventually concatenate the six datasets used here into one large dataset. However, we will keep them separate for now, as some of the steps to follow require individual datasets to be calibrated separately (namely, the sky/Tsys calibration and baseline subtraction). We therefore start by defining an array "basename" that includes the names of the six files in chronological order. This will simplify the following steps by allowing us to loop through the files using a simple for-loop in python. Remember that if you log out of CASA, you will have to re-issue this command. We will remind you of this in the relevant sections by repeating the command at the start.
# In CASA
basename=['uid___A002_X60b415_X39a','uid___A002_X60b415_X6f7','uid___A002_X6218fb_X264', 'uid___A002_X6218fb_X425','uid___A002_X6321c5_X3a7','uid___A002_X6321c5_X5ca']
Creating the MS and ASAP dataset
The raw data have been provided to you in the ASDM format. ASDM stands for ALMA Science Data Model. It is the native format of the data produced by the observatory. Before we can proceed to the calibration, we will need to convert those data to the CASA MS format. This is done simply with the task importasdm. For example:
In CASA
importasdm(vis = 'uid___A002_X60b415_X39a')
Note: importasdm has an option singledish, which you may be tempted to use. It works, but it has some limitations (which will be removed in the future), so for now, we recommend not using it.
The usual first step is then to get some basic information about the data. We do this using the task listobs, which will output a detailed summary of each dataset supplied.
# In CASA
for name in basename:
listobs(vis=name+'.ms')
Note that after cutting and pasting a for-loop you often have to press return several times to execute. The output will be sent to the CASA logger. You will have to scroll up to see the individual output for each of the six datasets. Here is an example of the most relevant output for the fifth file in the list.
Observation: ALMA
Data records: 16588 Total integration time = 1302.91 seconds
Observed from 12-May-2013/02:22:16.9 to 12-May-2013/02:43:59.8 (UTC)
ObservationID = 0 ArrayID = 0
Date Timerange (UTC) Scan FldId FieldName nRows nUnflRows SpwIds Average Interval(s) ScanIntent
12-May-2013/02:22:16.3 - 02:22:42.8 1 0 M100 900 0.00 [0, 1, 2, 3, 4, 5, 6, 7, 8] [1.15, 0.48, 0.48, 0.48, 0.48, 0.48, 0.48, 0.48, 0.48] CALIBRATE_ATMOSPHE
RE#ON_SOURCE,CALIBRATE_WVR#ON_SOURCE
02:23:24.3 - 02:25:11.4 2 0 M100 1524 0.00 [0, 9, 10, 11, 12, 13, 14, 15, 16] [1.15, 1.01, 1.01, 1.01, 1.01, 1.01, 1.01, 1.01, 1.01] OBSERVE_TAR
GET#OFF_SOURCE,CALIBRATE_WVR#OFF_SOURCE
02:25:25.2 - 02:27:11.2 3 0 M100 1522 0.00 [0, 9, 10, 11, 12, 13, 14, 15, 16] [1.15, 1.01, 1.01, 1.01, 1.01, 1.01, 1.01, 1.01, 1.01] OBSERVE_TAR
GET#OFF_SOURCE,CALIBRATE_WVR#OFF_SOURCE
02:27:25.1 - 02:29:11.0 4 0 M100 1522 0.00 [0, 9, 10, 11, 12, 13, 14, 15, 16] [1.15, 1.01, 1.01, 1.01, 1.01, 1.01, 1.01, 1.01, 1.01] OBSERVE_TAR
GET#OFF_SOURCE,CALIBRATE_WVR#OFF_SOURCE
02:29:24.9 - 02:31:12.0 5 0 M100 1524 0.00 [0, 9, 10, 11, 12, 13, 14, 15, 16] [1.15, 1.01, 1.01, 1.01, 1.01, 1.01, 1.01, 1.01, 1.01] OBSERVE_TAR
GET#OFF_SOURCE,CALIBRATE_WVR#OFF_SOURCE
02:31:25.8 - 02:33:11.8 6 0 M100 1522 0.00 [0, 9, 10, 11, 12, 13, 14, 15, 16] [1.15, 1.01, 1.01, 1.01, 1.01, 1.01, 1.01, 1.01, 1.01] OBSERVE_TAR
GET#OFF_SOURCE,CALIBRATE_WVR#OFF_SOURCE
02:34:03.6 - 02:34:29.0 7 0 M100 900 0.00 [0, 1, 2, 3, 4, 5, 6, 7, 8] [1.15, 0.48, 0.48, 0.48, 0.48, 0.48, 0.48, 0.48, 0.48] CALIBRATE_ATMOSPHE
RE#ON_SOURCE,CALIBRATE_WVR#ON_SOURCE
02:34:42.8 - 02:36:30.0 8 0 M100 1524 0.00 [0, 9, 10, 11, 12, 13, 14, 15, 16] [1.15, 1.01, 1.01, 1.01, 1.01, 1.01, 1.01, 1.01, 1.01] OBSERVE_TAR
GET#OFF_SOURCE,CALIBRATE_WVR#OFF_SOURCE
02:36:43.8 - 02:38:29.8 9 0 M100 1524 0.00 [0, 9, 10, 11, 12, 13, 14, 15, 16] [1.15, 1.01, 1.01, 1.01, 1.01, 1.01, 1.01, 1.01, 1.01] OBSERVE_TAR
GET#OFF_SOURCE,CALIBRATE_WVR#OFF_SOURCE
02:38:43.6 - 02:40:30.7 10 0 M100 1524 0.00 [0, 9, 10, 11, 12, 13, 14, 15, 16] [1.15, 1.01, 1.01, 1.01, 1.01, 1.01, 1.01, 1.01, 1.01] OBSERVE_TAR
GET#OFF_SOURCE,CALIBRATE_WVR#OFF_SOURCE
02:40:44.5 - 02:42:30.5 11 0 M100 1524 0.00 [0, 9, 10, 11, 12, 13, 14, 15, 16] [1.15, 1.01, 1.01, 1.01, 1.01, 1.01, 1.01, 1.01, 1.01] OBSERVE_TAR
GET#OFF_SOURCE,CALIBRATE_WVR#OFF_SOURCE
02:42:44.4 - 02:44:00.4 12 0 M100 1078 0.00 [0, 9, 10, 11, 12, 13, 14, 15, 16] [1.15, 1.01, 1.01, 1.01, 1.01, 1.01, 1.01, 1.01, 1.01] OBSERVE_TAR
GET#OFF_SOURCE,CALIBRATE_WVR#OFF_SOURCE
(nRows = Total number of rows per scan)
Fields: 1
ID Code Name RA Decl Epoch SrcId nRows nUnflRows
0 none M100 12:22:54.899040 +15.49.20.57200 J2000 0 16588 0.00
Spectral Windows: (17 unique spectral windows and 2 unique polarization setups)
SpwID Name #Chans Frame Ch1(MHz) ChanWid(kHz) TotBW(kHz) BBC Num Corrs
0 WVR#NOMINAL 4 TOPO 184550.000 1500000.000 7500000.0 0 XX
1 ALMA_RB_03#BB_1#SW-01#FULL_RES 128 TOPO 101942.187 -15625.000 2000000.0 1 XX YY
2 ALMA_RB_03#BB_1#SW-01#CH_AVG 1 TOPO 100926.562 1781250.000 1781250.0 1 XX YY
3 ALMA_RB_03#BB_2#SW-01#FULL_RES 128 TOPO 103757.337 -15625.000 2000000.0 2 XX YY
4 ALMA_RB_03#BB_2#SW-01#CH_AVG 1 TOPO 102741.712 1781250.000 1781250.0 2 XX YY
5 ALMA_RB_03#BB_3#SW-01#FULL_RES 128 TOPO 111814.962 15625.000 2000000.0 3 XX YY
6 ALMA_RB_03#BB_3#SW-01#CH_AVG 1 TOPO 112783.712 1781250.000 1781250.0 3 XX YY
7 ALMA_RB_03#BB_4#SW-01#FULL_RES 128 TOPO 113689.962 15625.000 2000000.0 4 XX YY
8 ALMA_RB_03#BB_4#SW-01#CH_AVG 1 TOPO 114658.712 1781250.000 1781250.0 4 XX YY
9 ALMA_RB_03#BB_1#SW-01#FULL_RES 4080 TOPO 101945.850 -488.281 1992187.5 1 XX YY
10 ALMA_RB_03#BB_1#SW-01#CH_AVG 1 TOPO 100949.756 1992187.500 1992187.5 1 XX YY
11 ALMA_RB_03#BB_2#SW-01#FULL_RES 4080 TOPO 103761.000 -488.281 1992187.5 2 XX YY
12 ALMA_RB_03#BB_2#SW-01#CH_AVG 1 TOPO 102764.906 1992187.500 1992187.5 2 XX YY
13 ALMA_RB_03#BB_3#SW-01#FULL_RES 4080 TOPO 111811.300 488.281 1992187.5 3 XX YY
14 ALMA_RB_03#BB_3#SW-01#CH_AVG 1 TOPO 112806.906 1992187.500 1992187.5 3 XX YY
15 ALMA_RB_03#BB_4#SW-01#FULL_RES 4080 TOPO 113686.300 488.281 1992187.5 4 XX YY
16 ALMA_RB_03#BB_4#SW-01#CH_AVG 1 TOPO 114681.906 1992187.500 1992187.5 4 XX YY
Sources: 19
ID Name SpwId RestFreq(MHz) SysVel(km/s)
0 M100 0 - -
0 M100 17 - -
0 M100 18 - -
0 M100 1 - -
0 M100 2 - -
0 M100 3 - -
0 M100 4 - -
0 M100 5 - -
0 M100 6 - -
0 M100 7 - -
0 M100 8 - -
0 M100 9 100950 0
0 M100 10 100950 0
0 M100 11 102794.1 0
0 M100 12 102794.1 0
0 M100 13 112794.1 0
0 M100 14 112794.1 0
0 M100 15 114669.1 0
0 M100 16 114669.1 0
Antennas: 2:
ID Name Station Diam. Long. Lat. Offset from array center (m) ITRF Geocentric coordinates (m)
East North Elevation x y z
0 PM01 T704 12.0 m -067.45.16.2 -22.53.22.1 42.8992 -520.1885 22.2159 2225113.044955 -5440122.820877 -2481517.728410
1 PM04 T703 12.0 m -067.45.16.2 -22.53.23.9 42.8809 -575.6904 21.7744 2225104.701420 -5440102.470120 -2481568.688227
We are now going to go through each of them with a bit more explanations. We will take the example of uid___A002_X6218fb_X264.ms. Before starting, you need to know that most of the tasks that we will use are part of the ASAP package, which was incorporated into CASA. The ASAP package is using a different data format, so from a global point of view, what we are going to do is, first convert the MS to the ASAP format, then run the necessary calibration tasks, then convert the data back to the MS format. Another important difference with interferometric data reduction is that the calibration is performed directly on the dataset, we will not produce calibration tables and apply them at the end. An effort is on-going to update the SD routines so that this is done, that should be available soon, but until then, please remember that all SD calibration operations apply to the data directly, so you may want to always create a new dataset each time, so that you do not have to start all over again.
We will then convert the data to the ASAP format. We will use the routine sd.splitant for that. An outprefix has to be specified because the routine will create a dataset per antenna, taking the outprefix, and appending the antenna name and the format extension ('.asap').
# In CASA
sd.splitant(filename = 'uid___A002_X6218fb_X264.ms',
outprefix = 'uid___A002_X6218fb_X264.ms',
overwrite = True,
freq_tolsr = True)
The reason for the option freq_tolsr=True is to convert the frequencies in the MS from TOPO to LSRK. This is an important setting, as due to current limitations, it is the only stage where this conversion can be done. Using freq_tolsr=False would keep the frame of the frequencies from the MS, and would force you to image the final cube in TOPO velocities. Another important aspect is for the baselining/line finding: if your observations have been obtained at different epochs, it may be easier to work on spectra with LSRK frequencies because then the line emissions of all spectra can be overlapped.
Calibration
About the calibration itself: the two main steps are
1. the calibration of the spectra into K, by applying the Tsys calibration and removing the signal from the OFF position 2. removing the baselines (i.e. subtracting the background emission, to keep only the line emission.)
We have now two ASAP datasets, one for PM03 and one for PM04. As usual, we will first obtain information about the content of the datasets, using sdlist (equivalent of listobs).
# In CASA
sdlist(infile = 'uid___A002_X6218fb_X264.ms.PM03.asap',
outfile = 'uid___A002_X6218fb_X264.ms.PM03.asap.sdlist')
sdlist(infile = 'uid___A002_X6218fb_X264.ms.PM04.asap',
outfile = 'uid___A002_X6218fb_X264.ms.PM04.asap.sdlist')
Tsys
Let's start by checking the Tsys solutions. We will use the gencal command used for interferometric data reduction. This will produce a CASA calibration table, which can be analysed similarly to previous guides M100 Band3 ACA 4.1#Tsys
# In CASA
gencal(vis = 'uid___A002_X6218fb_X264.ms',
caltable = 'uid___A002_X6218fb_X264.ms.tsys',
caltype = 'tsys')
plotbandpass(caltable='uid___A002_X6218fb_X264.ms.tsys', overlay='time',
xaxis='freq', yaxis='amp', subplot=22, buildpdf=False, interactive=False,
showatm=True,pwv='auto',chanrange='5~123',showfdm=True,
field='', figfile='uid___A002_X6218fb_X264.ms.tsys.plots.overlayTime/uid___A002_X6218fb_X264.ms.tsys')
This sequence loops over all of our files and plots Tsys as a function of time for channel. In the call to plotcal:
- subplot=22 parameter sets up a 2 x 2 panel grid.
- iteration tells plotcal to make a separate plot for each antenna.
For the spectral window association, we will use the routine tsysspwmap available directly in CASA.
# In CASA
from recipes.almahelpers import tsysspwmap
tsysmap = tsysspwmap(vis = 'uid___A002_X6218fb_X264.ms', tsystable = 'uid___A002_X6218fb_X264.ms.tsys')
tsysmap is a list of indices, with tsysmap[i] being the index of the Tsys spw to associate to the science spw of index i.
In the next step, we will do the actual preparation of the Tsys (i.e. re-sampling and re-indexing of the Tsys spw), using the library filltsys, also available directly in CASA.
# In CASA
import filltsys
for i in [9, 11, 13, 15]:
filltsys.fillTsys('uid___A002_X6218fb_X264.ms.PM03.asap',
specif = i,
tsysif = tsysmap[i],
mode = 'linear',
extrap = True)
for i in [9, 11, 13, 15]:
filltsys.fillTsys('uid___A002_X6218fb_X264.ms.PM04.asap',
specif = i,
tsysif = tsysmap[i],
mode = 'linear',
extrap = True)
This will loop through each science spw, re-sample the Tsys solutions associated to spw of index tsysmap[i] using a linear interpolation, allowing for extrapolation in case the science and Tsys spw do not perfectly match, and re-index it to spw of index i. You will note that the SD tasks use the word 'if': this is equivalent to 'spw id' in interferometry tasks.
We have now prepared the Tsys solutions for application. We have not applied them yet. We will first do some a-priori flagging, of the edge channels. This is similar to flagging edge channels in TDM interferometry datasets. In FDM interferometry datasets, these edge channels are not written out, so there is no need to flag anything. Here, we have an FDM SD dataset, produced with the second correlator (the ACA correlator), which does not excise them (yet).
Each spw covers the full baseband width (2GHz), so we will flag 120 channels on each side so as to keep only the center 3840 channels, similarly to FDM interferometry datasets.
# In CASA
sdflag(infile = 'uid___A002_X6218fb_X264.ms.PM03.asap',
specunit = 'channel',
iflist = [9, 11, 13, 15],
maskflag = [[0, 119], [3960, 4079]],
overwrite = True)
sdflag(infile = 'uid___A002_X6218fb_X264.ms.PM04.asap',
specunit = 'channel',
iflist = [9, 11, 13, 15],
maskflag = [[0, 119], [3960, 4079]],
overwrite = True)
Now is the step where we will actually calibrate the data. This is done with the task sdcal. This is a straight-forward operation, so no opportunity for tweaking. The calibration mode (calmode) shall be 'ps', as in Position Switching. The task will find automatically the OFF-position measurements. We need to specify the science spws (iflist) because they are the only ones that can be calibrated. We specify a different output file (outfile) so as to be able to come back to the original dataset if necessary.
# In CASA
sdcal(infile = 'uid___A002_X6218fb_X264.ms.PM03.asap',
calmode = 'ps',
iflist = [9, 11, 13, 15],
scanaverage = False,
timeaverage = False,
polaverage = False,
outfile = 'uid___A002_X6218fb_X264.ms.PM03.asap.cal',
overwrite = True)
sdcal(infile = 'uid___A002_X6218fb_X264.ms.PM04.asap',
calmode = 'ps',
iflist = [9, 11, 13, 15],
scanaverage = False,
timeaverage = False,
polaverage = False,
outfile = 'uid___A002_X6218fb_X264.ms.PM04.asap.cal',
overwrite = True)
Before proceeding to the next step, we can plot the calibrated spectra, using the sdplot task. The commands below will plot one spectrum per scan, spw and polarization. It is difficult to describe what a good spectrum looks like at this stage. What you should look for is outliers (i.e. spectra that differ from most, whether in shape or level).
# In CASA
for i in [9, 11, 13, 15]:
sdplot(infile='uid___A002_X6218fb_X264.ms.PM03.asap.cal',
iflist=[i], plottype='spectra', specunit='channel', scanaverage=True, stack='pol', panel='scan',
outfile='uid___A002_X6218fb_X264.ms.PM03.asap.cal.plots/uid___A002_X6218fb_X264.ms.PM03.asap.cal.spectra.spw'+str(i)+'.png',
overwrite = True)
for i in [9, 11, 13, 15]:
sdplot(infile='uid___A002_X6218fb_X264.ms.PM04.asap.cal',
iflist=[i], plottype='spectra', specunit='channel', scanaverage=True, stack='pol', panel='scan',
outfile='uid___A002_X6218fb_X264.ms.PM04.asap.cal.plots/uid___A002_X6218fb_X264.ms.PM04.asap.cal.spectra.spw'+str(i)+'.png',
overwrite = True)
Convert each ASAP dataset to an MS
In CASA
sdsave(infile = 'uid___A002_X6218fb_X264.ms.PM03.asap.cal.bl',
outfile = 'uid___A002_X6218fb_X264.ms.PM03.asap.cal.bl.ms',
outform = 'MS2')
sdsave(infile = 'uid___A002_X6218fb_X264.ms.PM04.asap.cal.bl',
outfile = 'uid___A002_X6218fb_X264.ms.PM04.asap.cal.bl.ms',
outform = 'MS2')
Combine all executions to one MS
Concatenate all of the calibrated measurement sets into one for imaging. The CASA task "concat" will do this.
In CASA
os.system('rm -rf concat_m100.ms')
concat(vis='uid___A002_X60b415_X39a.ms.cal.split', 'uid___A002_X60b415_X6f7.ms.cal.split', 'uid___A002_X6218fb_X264.ms.cal.split', 'uid___A002_X6218fb_X425.ms.cal.split', 'uid___A002_X6321c5_X3a7.ms.cal.split', 'uid___A002_X6321c5_X5ca.ms.cal.split'],
concatvis='concat_m100.ms',
freqtol='10MHz')
The individual calibrated MSs have slightly different observing frequencies, although the rest frequencies are the same. The freqtol parameter sets the tolerance for considering whether the different spectral windows from the input datasets should be output as the same spectral window ID.
Image the Total Power Data
Run listobs on the total power data to see what spw contains the CO
In CASA
os.system('rm -rf concat_m100.ms.listobs')
listobs(vis='concat_m100.ms',listfile='concat_m100.ms.listobs')
Spectral window SPWID=3 contains the 115.27 GHz line, so we image this window. The task "sdimaging" will do this.
In CASA
os.system('rm -rf TP_CO_cube')
sdimaging(infile='concat_m100.ms',
field=0,spw=3,
specunit='km/s',restfreq='115.271204GHz',
dochannelmap=True,
nchan=70,start=1400,step=5,
gridfunction='gjinc',imsize=[50,50],
cell=['10arcsec','10arcsec'],
outfile='TP_CO_cube')
The restfreq parameter must be specified when using "km/s" as the units, as in this case. Start and step parameters are specified in units that the user chooses for specunit. The numbers here are chosen so that the resulting image has the same number of channels, velocity range and channel width as the 7m and 12m array images. The gridfunction is the weighting function that is used to grid the observed flux to individual pixels in the image. "SF" is a spheroidal function, which minimizes aliasing effects. "BOX" is a pillbox function, which defaults to a kernel box size of 1 pixel. The "PB" (primary beam) assumes an Airy disk, corresponding to an antenna with 10.7m diameter, the effective diameter of an ALMA 12m antenna. The "GAUSS" is a gaussian, and its size can be defined by additional subparameters (truncate and gwidth). "GJINC" is a gaussian convolved with the Bessel function, and can minimize the broadening of the effective beam. Any of the functions which require the obseving frequency for determining the beam size will read the frequency from the dataset, and the user can use the default.
The cell size should be chosen so that it is about 1/3 to 1/4 of the FWHM of the effective beam.
