TWHydraBand7 Calibration for CASA 3.3
- This script assumes that you have downloaded RawDataAndTablesForReduction from TWHydraBand7#Getting_the_Data
- Details of the ALMA data are provided at TWHydraBand7
Unpack the Data
Once you have downloaded RawDataAndTablesForReduction.tar.gz, in a terminal unpack it with
>tar -zxvf RawDataAndTablesForReduction.tar.gz
It is a good idea to rename the directory in case you download another SV dataset
>mv RawDataAndTablesForReduction TWHyaBand7_RawDataAndTablesForReduction
cd TWHyaBand7_RawDataAndTablesForReduction
start CASA
>casapy
Observing Log and Priors
Below is some information that may not be obvious from the data
Like at most telescopes, the ALMA operators and astronomers at the site log any significant issues that were noticed from the online system during the observation. For these Band 7 observations of TW Hya on April 22, 2011, the only issue of note was that antenna DV04 could not observe at Band 7 and will need to be flagged.
First lets examine the basic properties of the data using listobs. Among other information we will be able to deduce what the source ids are and what each source was used for:
3c279 is the bandpass calibrator (field id=0) Titan is the absolute flux calibrator (field id=1) J1147-382 is the secondary phase calibrator (field id=3) J1037-295 is the primary phase calibrator field id=4)
# In CASA
data=['X3c1.ms','X5d8.ms','X7ef.ms']
for vis in data:
listobs(vis=vis,verbose=T)
Below we copy select parts of the listobs output reported in the CASA logger for reference.
X3c1.ms: Data records: 278352 Total integration time = 6109.06 seconds Observed from 22-Apr-2011/00:01:52.9 to 22-Apr-2011/01:43:42.0 (UTC) X5d8.ms: Data records: 278406 Total integration time = 5995.01 seconds Observed from 22-Apr-2011/01:48:05.8 to 22-Apr-2011/03:28:00.8 (UTC) X7ef.ms: Data records: 255717 Total integration time = 5407.49 seconds Observed from 22-Apr-2011/03:30:39.7 to 22-Apr-2011/05:00:47.2 (UTC) ID Code Name RA Decl Epoch SrcId 0 none 3c279 12:56:11.1666 -05.47.21.5247 J2000 0 1 none Titan 00:00:00.0000 +00.00.00.0000 J2000 1 2 none TW Hya 11:01:51.8450 -34.42.17.1609 J2000 2 3 none J1147-382=Q* 11:47:01.3815 -38.12.11.1179 J2000 3 4 none J1037-295=Q* 10:37:16.0899 -29.34.02.9888 J2000 4 SpwID Chans Frame Ch1(MHz) ChanWid(kHz)TotBW(kHz) Ref(MHz) Corrs 0 4 TOPO 184550 1500000 7500000 183300 I 1 128 TOPO 355740.062 15625 2000000 355732.25 XX YY 2 1 TOPO 356716.625 1796875 1796875 355732.25 XX YY 3 128 TOPO 356507.813 15625 2000000 356500 XX YY 4 1 TOPO 357484.375 1796875 1796875 356500 XX YY 5 128 TOPO 346792.187 15625 2000000 346800 XX YY 6 1 TOPO 345784.375 1796875 1796875 346800 XX YY 7 128 TOPO 345182.438 15625 2000000 345190.25 XX YY 8 1 TOPO 344174.625 1796875 1796875 345190.25 XX YY 9 128 TOPO 344386.763 15625 2000000 344394.575 XX YY 10 1 TOPO 343378.95 1796875 1796875 344394.575 XX YY 11 128 TOPO 346324.263 15625 2000000 346332.075 XX YY 12 1 TOPO 345316.45 1796875 1796875 346332.075 XX YY 13 128 TOPO 354402.388 15625 2000000 354394.575 XX YY 14 1 TOPO 355378.95 1796875 1796875 354394.575 XX YY 15 128 TOPO 356402.388 15625 2000000 356394.575 XX YY 16 1 TOPO 357378.95 1796875 1796875 356394.575 XX YY 17 3840 TOPO 356497.936 122.070312 468750 355732.25 XX YY 18 1 TOPO 356732.189 468750 468750 355732.25 XX YY 19 3840 TOPO 357734.314 122.070312 468750 356500 XX YY 20 1 TOPO 357499.939 468750 468750 356500 XX YY 21 3840 TOPO 346034.314 122.070312 468750 346800 XX YY 22 1 TOPO 345799.939 468750 468750 346800 XX YY 23 3840 TOPO 343955.936 122.070312 468750 345190.25 XX YY 24 1 TOPO 344190.189 468750 468750 345190.25 XX YY ID= 0-3: 'DV04'='J505', 'DV06'='T704', 'DV07'='J510', 'DV08'='T703', ID= 4-7: 'DV09'='N602', 'DV10'='N606', 'PM01'='T702', 'PM02'='T701', ID= 8-8: 'PM03'='J504'
Things to Notice from above:
- Titan appears to have a position of 0,0. We will fix this later in the tutorial (it is a temporary data capture issue in the ALMA system, Titan was properly tracked during the observation).
- There are a lot of spectral windows! Don't be alarmed. In ALMA data,
- spw=0 is always reserved for the water vapor radiometry data.
- spw=1, 3, 5, 7 are the wideband (TDM) correlator mode with 128 channels per spw default settings for doing pointing
- currently Tsys data can only be taken in the wide (TDM) correlator mode with 128 channels per spw, these correspond to spws=9, 11, 13, 15 above. The frequencies of these TDM spw were carefully chosen so the the narrowband FDM spws fall within them.
- The 0.5 GHz, high spectral resolution spws that contain all the real source data (science target and calibrators) are 17, 19, 21, and 23
- For quicklook purposes, the ALMA online system produces a channel averaged spw for each channelized spw, these correspond to spw=2,4,6,8,10,12,14,16,18,20,22,24 and shouldn't be used for anything in post-processing as bandpass calibration etc has not been applied.
Next scroll up in the logger itself after making your logger window as wide as possible. In the portion of the listobs that shows the timerange, the final column gives the scan intent. It is important to have a look at this because you can flag and do other operations in CASA using the scan intent as a selection option. Also these intents will be used in the future for pipeline processing. You will usually want to use them to flag the CALIBRATE_POINTING* and CALIBRATE_ATMOS* (the latter correspond to the hot and cold load measurements used to calculate Tsys) data. This is essential for science or calibration data taken only in TDM (wide bandwidth) mode, i.e. the same mode that the pointing and Tsys data are taken in. For the current data however, that calibration data is in different spws from the useful FDM (narrow bandwidth) data, so we'll separate them by splitting them off below.
Examine Tsys and WVR Tables, Apply, and Split
A number of tables were included in the RawDataAndTablesForReduction directory, including ones of the form *.tsys, *.tsys.fdm, and wvr*cal. The *.tsys tables are the Tsys values derived from the TDM (wide bandwidth) correlator mode data. They were used to derive Tsys for their corresponding FDM (narrow bandwidth) spws: *.tsys.fdm. The wvr*cal contain the phase corrections derived from the Water Vapor Radiometers. These tables have been provided to you because they are currently generated from scripts rather than inside CASA, a situation that will change in Early Science.
It is very important that Tsys measurements as close in time and elevation to a particular source are applied. To save time, Tsys measurements are sometimes skipped, so the first step is to establish which sources have their own Tsys measurements and which may not. To do this we will make a plot of X3c1.tsys for each source; if the source has Tsys measurements you will see them (the two polarizations are different colors). We select only one channel to speed plotting.
# In CASA
plotcal(caltable='X3c1.tsys',xaxis='time',yaxis='amp',field='0',antenna='1~8',
iteration='antenna',subplot=421,poln='',spw='1:50~50')
# In CASA
plotcal(caltable='X3c1.tsys',xaxis='time',yaxis='amp',field='1',antenna='1~8',
iteration='antenna',subplot=421,poln='',spw='1:50~50')
# In CASA
plotcal(caltable='X3c1.tsys',xaxis='time',yaxis='amp',field='2',antenna='1~8',
iteration='antenna',subplot=421,poln='',spw='1:50~50')
# In CASA
plotcal(caltable='X3c1.tsys',xaxis='time',yaxis='amp',field='3',antenna='1~8',
iteration='antenna',subplot=421,poln='',spw='1:50~50')
For this one you will get a message like:
2011-05-24 19:51:55 SEVERE Exception Reported: Combined selection selects nothing. .*** Error *** Combined selection selects nothing.
Indicating that there are no Tsys measurements for this source.
# In CASA
plotcal(caltable='X3c1.tsys',xaxis='time',yaxis='amp',field='4',antenna='1~8',
iteration='antenna',subplot=421,poln='',spw='1:50~50')
So except for Field 3 (the secondary phase calibrator ), the other sources including TW Hya itself have their own Tsys measurements. Since the secondary phase calibrator is close on the sky to the primary phase calibrator and observed near in time, we will transfer the Tsys from the primary to the secondary phase calibrator. For sources with Tsys measurements, use applycal field and gainfield parameters to apply solutions only to themselves. While applying Tsys, we will also apply the wvr calibration tables; wvr data is usually always present for all sources.
Below we set some definitions so we can loop over the the correct parameters.
# In CASA
data=['X3c1.ms','X5d8.ms','X7ef.ms']# This source had no Tsys measurements associated with it
fdm='17,19,21,23'
sources=['0','1','2','4']
nocal='3'
# In CASA
for vis in data:
for field in sources:
applycal(vis=vis,spw=fdm,field=field,gainfield=field,
gaintable=['%s.tsys.fdm'%(vis.split('.')[0]),'wvr_%s.cal'%(vis.split('.')[0])],
interp=['nearest','nearest'],flagbackup=F)
For sources without Tsys (i.e. field=3), couple them with best source with Tsys (i.e. field=4).
# In CASA
for vis in data:
applycal(vis=vis,spw=fdm,field=nocal,
gaintable=['%s.tsys.fdm'%(vis.split('.')[0]),'wvr_%s.cal'%(vis.split('.')[0])],
gainfield=['4',nocal],interp=['nearest','nearest'],flagbackup=F)
Check Tsys application with plotms. First we make plots with ydatacolumn='data to look at the raw data. For simplicity we'll just look at the two strongest sources 3C279 and Titan.
# In CASA
plotms(vis='X5d8.ms',spw='17,19,21,23',xaxis='frequency',yaxis='amp',field='0,1',
antenna='*&*',avgtime='1e8',avgscan=T,coloraxis='field',iteraxis='spw',
xselfscale=T,ydatacolumn='data')
Note: the antenna='*&*' is specifying that we only want to see the cross-correlation data (i.e. not the autocorrelations).
Use the green arrows on the plotms display to cycle through the 4 FDM spws. Now change Data Column to corrected on the "Axes" tab (this is the same as ydatacolumn in the task) and hit plot button on the gui. Flip through again.
Things to notice:
- The amplitude scale has changed due to application of Tsys
- spw=17 is clean and flat
- spw=19 shows a strong dip around 358.1 GHz -- this is an atmospheric absorption feature. In the future high spectral resolution Tsys measurements may help take out some of the absorption, until then we will need to flag this part of the spectrum.
- spw=21 Titan shows both a narrow and a broad spectral feature. The broad feature is CO(3-2) and it completely fills the 0.5 GHz band so there are no line-free channels! Titan must be handled carefully when it is used as a calibrator because it contains strong line emission. Indeed, we will not be able to use the absolute flux transfer from this spw for other targets. When possible Titan should be avoided for absolute flux calibration. For these observations, no other suitable absolute flux calibrator was available.
- spw=23 shows another strong line from Titan.
Beyond these things to notice, the application of Tsys and wvr seem to have gone well. We can now split off only the narrow band FDM spectral windows for further processing.
# In CASA
for vis in data:
split(vis=vis,outputvis='%s_wvrtsys.ms'%(vis.split('.')[0]),datacolumn='corrected',spw=fdm)
Note: After split, the spectral windows will be renumbered to 0, 1, 2, 3 corresponding to the old spw=17, 19, 21, 23, see listobs excerpt below.
SpwID #Chans Frame Ch1(MHz) ChanWid(kHz)TotBW(kHz) Ref(MHz) Corrs 0 3840 TOPO 356497.936 122.070312 468750 355732.25 XX YY 1 3840 TOPO 357734.314 122.070312 468750 356500 XX YY 2 3840 TOPO 346034.314 122.070312 468750 346800 XX YY 3 3840 TOPO 343955.936 122.070312 468750 345190.25 XX YY
Initial Inspection and Flagging
From here we begin to operate on the *wvrtsys.ms files.
ALMA data contains both the cross correlation and autocorrelation data. Presently nothing more can be done with the autocorrelation data so we flag it. Additionally, for smaller configurations of the array, and northerly sources one antenna can shadow another, blocking its view. This data also needs to be flagged. Finally, the observing log told us DV04 was not behaving properly at Band 7 and shouldn't be used so we flag that as well.
# In CASA
data=['X3c1_wvrtsys.ms','X5d8_wvrtsys.ms','X7ef_wvrtsys.ms']
for vis in data:
flagdata(vis=vis,autocorr=True, flagbackup=T)
flagdata(vis=vis,mode='shadow',diameter=12.0, flagbackup=F)
flagdata(vis=vis,antenna='DV04', flagbackup=F)
Next we inspect the data with plotms for additional issues that need flagging, first looking at amplitude vs. time for each of the three datasets. Use the green arrows on the plotms gui to cycle through spws.
Note: when iteraxis is invoked in plotms, unzoom does not always work properly. If you have trouble with this, flip to next page and then back to current one. this usually fixes the problem.
# In CASA
plotms(vis='X3c1_wvrtsys.ms',spw='',xaxis='time',yaxis='amp',field='',avgchannel='3840',
coloraxis='field',iteraxis='spw')
# In CASA
plotms(vis='X5d8_wvrtsys.ms',spw='',xaxis='time',yaxis='amp',field='',avgchannel='3840',
coloraxis='field',iteraxis='spw')
# In CASA
plotms(vis='X7ef_wvrtsys.ms',spw='',xaxis='time',yaxis='amp',field='',avgchannel='3840',
coloraxis='field',iteraxis='spw')
Things to Notice:
- In spw=2 there are some low points in all three data files. Use the "Mark Regions" box (left of center on bottom row of icons) to draw small box on low points. Then click the magnifying glass icon to the right of center. The output will be sent to the logger window. Note: keep the marked region small or you can overwhelm the buffer. Locate suggests correlation YY on PM03 is the culprit.
To see this more clearly you can rerun plotms with coloraxis='field' switched to coloraxis='corr'
# In CASA
plotms(vis='X3c1_wvrtsys.ms',spw='2',xaxis='time',yaxis='amp',field='',avgchannel='3840',
coloraxis='corr')
To further check, you can exclude PM03 from the plot using antenna='!PM03'.
# In CASA
plotms(vis='X3c1_wvrtsys.ms',spw='2',xaxis='time',yaxis='amp',field='',avgchannel='3840',
coloraxis='corr',antenna='!PM03')
Having confirmed the problem antenna, we flag it below.
# In CASA
data=['X3c1_wvrtsys.ms','X5d8_wvrtsys.ms','X7ef_wvrtsys.ms']
for vis in data:
flagdata(vis=vis, flagbackup=T,
spw=['2'],
antenna=['PM03'],
correlation=['YY'])
Next we need to inspect the spectral properties of the data to look for issues.
First we look at phase as a function of frequency for a well behaved antenna like DV06. coloraxis='baselines' (note this parameter is called colorize in the gui, Display Tab) will show each baseline as a different color. Significant delay or baseline error will show up as phase slopes in these plots. Again you will want to page through the spectral windows.
# In CASA
plotms(vis='X3c1_wvrtsys.ms',spw='',xaxis='frequency',yaxis='phase',field='0',antenna='DV06',
avgtime='1e8',avgscan=T,coloraxis='baseline',iteraxis='spw',xselfscale=T)
# In CASA
plotms(vis='X5d8_wvrtsys.ms',spw='',xaxis='frequency',yaxis='phase',field='0',antenna='DV06',
avgtime='1e8',avgscan=T,coloraxis='baseline',iteraxis='spw',xselfscale=T)
# In CASA
plotms(vis='X7ef_wvrtsys.ms',spw='',xaxis='frequency',yaxis='phase',field='0',antenna='DV06',
avgtime='1e8',avgscan=T,coloraxis='baseline',iteraxis='spw',xselfscale=T)
Overall these plots look good with no large delays or baseline errors apparent.
With such high spectral resolution data, one often sees occasional "birdies" in the data. These are typically very narrow weak spectral features that are internally generated in the system. Over time as we track down issues the number of birdies should decrease, but you should always check. This is easiest to do by looking at amplitude vs. frequency for sources observed over a long time range (for sensitivity). It is helpful to look at more than one source to verify the nature of the emission. Below we look at both calibrators (brown and green) and TW Hya (orange). If you have enough sensitivity birdies should be seen on all sources - however, be careful not to mistake real line emission for a birdie even on a calibrator they can be present as we've seen for Titan.
# In CASA
plotms(vis='X3c1_wvrtsys.ms',spw='',xaxis='channel',yaxis='amp',field='2,3,4',
avgtime='1e8',avgscan=T,coloraxis='field',iteraxis='spw',xselfscale=T,yselfscale=T)
Do a quick check that other two datasets are the same (they are usually similar if the source of internal resonance hasn't changed in the meantime).
# In CASA
plotms(vis='X5d8_wvrtsys.ms',spw='',xaxis='channel',yaxis='amp',field='2,3,4',
avgtime='1e8',avgscan=T,coloraxis='field',iteraxis='spw',xselfscale=T,yselfscale=T)
# In CASA
plotms(vis='X7ef_wvrtsys.ms',spw='',xaxis='channel',yaxis='amp',field='2,3,4',
avgtime='1e8',avgscan=T,coloraxis='field',iteraxis='spw',xselfscale=T,yselfscale=T)
Now you can zoom in on each region to get channels for flagging, or to save time just go ahead and run the flagging command below where this has been done for you.
# In CASA
data=['X3c1_wvrtsys.ms','X5d8_wvrtsys.ms','X7ef_wvrtsys.ms']
for vis in data:
flagdata(vis=vis, flagbackup=T,
spw=['0:1067~1068;1279~1280;2367~2368;3775~3776',
'1:1279~1280;2367~2368;3775~3776',
'2:1279~1280;3775~3776',
'3:831~832;1535~1536;2367~2368;3775~3776;3839~3839'])
Replot in plotms to check that flagging has been done as anticipated.
Concatenate the Data
Here we will concatenate the three datasets for convience. WARNING: it is not always safe to do this [explain].
# In CASA
data=['X3c1_wvrtsys.ms','X5d8_wvrtsys.ms','X7ef_wvrtsys.ms']
concat(vis=data,concatvis='Band7multi_april22.ms')
If you like you can run listobs on new combined dataset
# In CASA
listobs(vis='Band7multi_april22.ms',verbose=F)
Fix Titan Coordinates
Note that in the near future this step will be unnecessary.
It is temporarily necessary to correct the positions of ephemeris objects observed by ALMA, as they are currently erroneously written to the data with coordinates of 00:00:00.0000 +00.00.00.0000 J2000 (see listobs output in TWHydraBand7_Calibration#Observing_Log_and_Priors). For this data we need to correct Titan which we observed as the absolute flux calibrator. It is important to apply the fix after concatenation of the data, or Titan will be at slightly different positions in all three datasets.
initialize CASA script with
# In CASA
execfile(casadef.python_library_directory+'/recipes/fixplanets.py')
Now fix the position of Titan. This step takes a while because the model and corrected scratch columns must be created for the new dataset first (see logger message).
# In CASA
fixplanets('Band7multi_april22.ms', 'Titan', True)
Check the header information again to see that the information for Titan has changed.
# In CASA
listobs(vis='Band7multi_april22.ms',verbose=F)
ID Code Name RA Decl Epoch SrcId 0 none 3c279 12:56:11.1666 -05.47.21.5247 J2000 0 1 none Titan 12:49:26.5120 -02.22.10.8190 J2000 1 2 none TW Hya 11:01:51.9063 -34.42.17.0212 J2000 2 3 none J1147-382=Q* 11:47:01.3815 -38.12.11.1179 J2000 3 4 none J1037-295=Q* 10:37:16.0899 -29.34.02.9888 J2000 4
Flag Calibrator Spectral Features
In this section we will flag spectral features in the calibrators so they do not skew the calibration solutions. First make plots of 3C279 and Titan in channel space to get channels for flagging.
# In CASA
plotms(vis='Band7multi_april22.ms',spw='',xaxis='channel',yaxis='amp',field='0',
avgtime='1e8',coloraxis='field',iteraxis='spw')
# In CASA
plotms(vis='Band7multi_april22.ms',spw='',xaxis='channel',yaxis='amp',field='1',
avgtime='1e8',coloraxis='field',iteraxis='spw')
Explicitly back up flag tables so flagged regions can be restored later.
# In CASA
flagmanager(vis='Band7multi_april22.ms',mode='save',versionname='before_calspectralflags')
Flag absorption features for all sources.
# In CASA
flagdata(vis='Band7multi_april22.ms',
spw=['1:2000~3840','2:1945~1960'])
Now we need to flag emission. The very broad CO(3-2) in spw=2 of Titan will prevent us from using it as a calibrator for this spectral window since there are no line-free channels (we will transfer, phase, amplitude, and absolute flux from another Titan spw below). We do want to flag the narrower Titan line in spw=3.
# In CASA
flagmanager(vis='Band7multi_april22.ms',mode='save',versionname='before_emissionflags')
# In CASA
flagdata(vis='Band7multi_april22.ms',
field=['1'],
spw=['3:1000~3000'])
Check if you like that the right things have been flagged.
# In CASA
plotms(vis='Band7multi_april22.ms',spw='',xaxis='channel',yaxis='amp',field='0,1',
avgtime='1e8',coloraxis='field',iteraxis='spw')
Set Up the Flux Calibrator Model
Fill the model data column for Titan with a model for the flux density as a function of baseline. The model is a uniformly illuminated disk with the size implied by the JPL Horizons ephemeris.
# In CASA
setjy(vis='Band7multi_april22.ms',field='1',
standard='Butler-JPL-Horizons 2010',scalebychan=F)
Bandpass Calibration
First we need to check how the amplitude and phase of 3C279 behave as a function of time.
# In CASA
plotms(vis='Band7multi_april22.ms',spw='',xaxis='time',yaxis='amp',
field='0',avgchannel='3840',iteraxis='spw')
Page through the spws. This shows significant variations in the amplitude across the three scans.
# In CASA
plotms(vis='Band7multi_april22.ms',spw='0',xaxis='time',yaxis='phase',
antenna='DV06;!DV04',coloraxis='corr',
field='0',avgchannel='3840',iteraxis='baseline')
Page through the baselines to DV06. This shows (1) significant variation in the phase across the three scans, as well as (2) smaller but significant variations within a scan. The first of these suggests that we do not want to combine the scans for bandpass calibration, but instead derive an independent bandpass for each scan. Actually what is more important than a DC offset is if the shape of the amplitude or phase change from scan to scan. There isn't enough S/N to tell for sure so we will keep them separate. The second suggest that we need to correct the phase vs. time behavior of the bandpass calibrator (3C279; field=0) within each scan, before doing the bandpass. We want to chose a relative narrow range of channels near the center of the spws so that bandpass phase slopes (to be corrected with the bandpass calibration itself) do not decorrelate the signal. You need to avoid spectral regions that were completely flagged above.
# In CASA
gaincal(vis='Band7multi_april22.ms',caltable='bpphase.gcal',
field='0',spw='0~3:900~1100',refant='DV06',
calmode='p',solint='int',minsnr=2.0,minblperant=4)
Inspect the phase-only calibration solutions for the bandpass calibrator. Look for any noisy solutions. We do XX and YY separately to reduce confusion in the plot. It will help to make your plotcal window wide
# In CASA
plotcal(caltable='bpphase.gcal',xaxis='time',yaxis='phase',spw='',antenna='1~8',
iteration='antenna',subplot=421,plotrange=[0,0,-180,180],
figfile='bpphase.X.png',poln='X')
In the plotcal gui you can use the magnifying glass symbol to zoom in on one of the scans to see the phase variations. The colors are the 4 different spws.
# In CASA
plotcal(caltable='bpphase.gcal',xaxis='time',yaxis='phase',spw='',antenna='1~8',
iteration='antenna',subplot=421,plotrange=[0,0,-180,180],
figfile='bpphase.Y.png',poln='Y')
Next we apply this phase-only correction on the fly while calculating the bandpass solutions.
# In CASA
bandpass(vis='Band7multi_april22.ms',caltable='bandpass.bcal',
field='0',spw='',combine='',refant='DV06',
solint='inf',solnorm=T,minblperant=4,fillgaps=17,
gaintable=['bpphase.gcal'])
You will see some error messages, these are just related to the flagged channels. We set fillgaps=17 to interpolate across the smaller regions of flagged channels (mesospheric and birdies).
A few words about solint and combine:
The use of solint='inf' in bandpass will derive one bandpass solution for each 3C279 scan, unless combine='scan' (which is the default). Here we set combine= explicitly so that it does not combine the scans into one bandpass. Likewise, field (spw) boundaries are only crossed if combine='field' (combine='spw'), the latter two are not generally good ideas for bandpass solutions.
Now inspect the phase and amplitude solutions. You are looking for excursions from smooth fits, or particularly noisy solutions compared to the others. These plots take a while because there are so many channels per spw (3840!).
# In CASA
plotcal(caltable='bandpass.bcal',xaxis='freq',yaxis='amp',antenna='1~8',
iteration='antenna',subplot=421,spw='0',figfile='bandpass.ampspw0.png')
# In CASA
plotcal(caltable='bandpass.bcal',xaxis='freq',yaxis='phase',antenna='1~8',
iteration='antenna',subplot=421,spw='0',figfile='bandpass.phasespw0.png')
The phase on DV06 is very small because it was used as the reference antenna.
These plots should be repeated for each spw. Be sure to rename the plots file each time before running.
Gain Calibration
Now that we have a bandpass solution to apply we can solve for the antenna-based phase and amplitude gain calibration. Since the phase changes on a much shorter timescale than the amplitude, we will solve for them separately. In particular, if the phase changes significantly over a scan time, the amplitude would be decorrelated, if the un-corrected phase were averaged over this timescale. Note that we re-solve for the gain solutions of the bandpass calibrator, so we can derive new solutions that are corrected for the bandpass shape. Since the bandpass calibrator will not be used again, this is not strictly necessary, but is useful to check its calibrated flux density for example.
# In CASA
gaincal(vis='Band7multi_april22.ms',caltable='intphase.gcal',
field='0,1,3,4',spw='0~3:40~3800',refant='DV06',
calmode='p',solint='int',minsnr=2.0,minblperant=4,
gaintable=['bandpass.bcal'])
Here solint='int' coupled with calmode='p' will derive a single phase solution for each 10 second integration. Note that the bandpass table is applied on-the-fly before solving for the phase solutions, however the bandpass is NOT applied to the data permanently until applycal is run later on.
Although solint='int' (i.e. the integration time of 10 seconds) is the best choice to apply before for solving for the amplitude solutions, it is not a good idea to use this to apply to the target. This is because the phase-scatter within a scan can dominate the interpolation between calibrator scans. Instead, we also solve for the phase on the scan time, solint='inf' (but combine=' ', since we want one solution per scan) for application to the target later on. Unlike the bandpass task, for gaincal, the default of the combine parameter is combine=' '.
# In CASA
gaincal(vis='Band7multi_april22.ms',caltable='scanphase.gcal',
field='0,1,3,4',spw='0~3:40~3800',refant='DV06',
calmode='p',solint='inf',minsnr=2.0,minblperant=4,
gaintable=['bandpass.bcal'])
Alternatively, instead of making a separate phase solution for application to the target, one can also run smoothcal to smooth the solutions derived on the integration time.
Next we apply the bandpass and solint='int' phase-only calibration solutions on-the-fly to derive amplitude solutions. Here the use of solint='inf', but combine=' ' will result in one solution per scan interval.
# In CASA
gaincal(vis='Band7multi_april22.ms',caltable='amp.gcal',
field='0,1,3,4',spw='0~3:40~3800',refant='DV06',
calmode='ap',solint='inf',minsnr=2.0,minblperant=4,
gaintable=['bandpass.bcal','intphase.gcal'])
Now carefully inspect all these solutions looking for discrepant solutions. If you see any you will need to flag them and rerun the calibration tasks above.
We make plots in X and Y separately for clarity of the plots. The colors are showing the spectral windows.
# In CASA
plotcal(caltable='intphase.gcal',xaxis='time',yaxis='phase',antenna='1~8',
spw='',field='0,1,3,4',iteration='antenna',subplot=421,
plotrange=[0,0,-180,180],poln='X',figfile='intphase_X.png')
# In CASA
plotcal(caltable='intphase.gcal',xaxis='time',yaxis='phase',antenna='1~8',
spw='',field='0,1,3,4',iteration='antenna',subplot=421,
plotrange=[0,0,-180,180],poln='Y',figfile='intphase_Y.png')
# In CASA
plotcal(caltable='scanphase.gcal',xaxis='time',yaxis='phase',antenna='1~8',
spw='',field='0,1,3,4',iteration='antenna',subplot=421,
plotrange=[0,0,-180,180],poln='X',figfile='scanphase_X.png')
# In CASA
plotcal(caltable='scanphase.gcal',xaxis='time',yaxis='phase',antenna='1~8',
spw='',field='0,1,3,4',iteration='antenna',subplot=421,
plotrange=[0,0,-180,180],poln='Y',figfile='scanphase_Y.png')
Since we have taken out the phase as best we can by applying the solint='int' phase-only solution, this plot will give a good idea of the residual phase error. If you see scatter of more than a few degrees here, you should consider going back and looking for more data to flag, particularly bad timeranges etc.
# In CASA
plotcal(caltable='amp.gcal',xaxis='time',yaxis='phase',antenna='1~8',
spw='',field='0,1,3,4',plotrange=[0,0,-1,1],
iteration='antenna',subplot=421,figfile='amp_phase.png')
These are very small, as they should be. Now look at the amplitude solutions.
# In CASA
plotcal(caltable='amp.gcal',xaxis='time',yaxis='amp',antenna='1~8',
iteration='antenna',subplot=421,spw='',poln='X',
plotrange=[0,0,0.0,0.3],figfile='amp_X.png')
# In CASA
plotcal(caltable='amp.gcal',xaxis='time',yaxis='amp',antenna='1~8',
iteration='antenna',subplot=421,spw='',poln='Y',
plotrange=[0,0,0.0,0.3],figfile='amp_Y.png')
Next we use the flux calibrator (whose flux density was set in setjy above) to derive the flux of the other calibrators. Note that the flux table REPLACES the amp.gcal in terms of future application of the calibration to the data, i.e. the flux table contains both the amp.gcal and flux scaling. Unlike the gain calibration steps, this is not an incremental table.
# In CASA
fluxscale(vis='Band7multi_april22.ms',caltable='amp.gcal',
fluxtable='flux.cal',reference='1',refspwmap=[0,1,3,3])
Its a good idea to note down the derived fluxscale values for your records.
Found reference field(s): Titan Found transfer field(s): 3c279 J1147-382=QSO J1037-295=QSO Spw=2 will be referenced to spw=3 Flux density for 3c279 in SpW=0 is: 10.7324 +/- 0.180055 (SNR = 59.6062, nAnt= 8) Flux density for 3c279 in SpW=1 is: 10.953 +/- 0.163625 (SNR = 66.9394, nAnt= 8) Flux density for 3c279 in SpW=2 (ref SpW=3) is: 10.7283 +/- 0.191615 (SNR = 55.9888, nAnt= 7) Flux density for 3c279 in SpW=3 is: 10.1286 +/- 0.143758 (SNR = 70.4556, nAnt= 8) Flux density for J1147-382=QSO in SpW=0 is: 0.913886 +/- 0.0106354 (SNR = 85.9287, nAnt= 8) Flux density for J1147-382=QSO in SpW=1 is: 0.734573 +/- 0.00852182 (SNR = 86.1991, nAnt= 8) Flux density for J1147-382=QSO in SpW=2 (ref SpW=3) is: 0.917729 +/- 0.0101124 (SNR = 90.7532, nAnt= 7) Flux density for J1147-382=QSO in SpW=3 is: 1.12946 +/- 0.0130193 (SNR = 86.7525, nAnt= 8) Flux density for J1037-295=QSO in SpW=0 is: 0.890222 +/- 0.00878729 (SNR = 101.308, nAnt= 8) Flux density for J1037-295=QSO in SpW=1 is: 0.708348 +/- 0.00675153 (SNR = 104.917, nAnt= 8) Flux density for J1037-295=QSO in SpW=2 (ref SpW=3) is: 0.896367 +/- 0.00856783 (SNR = 104.62, nAnt= 7) Flux density for J1037-295=QSO in SpW=3 is: 1.11563 +/- 0.0102939 (SNR = 108.378, nAnt= 8)
Now look at the flux solutions.
# In CASA
plotcal(caltable='flux.cal',xaxis='time',yaxis='amp',antenna='1~8',
iteration='antenna',subplot=421,spw='',poln='X',
plotrange=[0,0,0.0,0.3],figfile='flux_X.png')
# In CASA
plotcal(caltable='flux.cal',xaxis='time',yaxis='amp',antenna='1~8',
iteration='antenna',subplot=421,spw='',poln='Y',
plotrange=[0,0,0.0,0.3],figfile='flux_Y.png')
Apply Calibration and Inspect
Now we can use the flag manager to restore the spectral flagging we did on the calibrators.
# In CASA
flagmanager(vis='Band7multi_april22.ms',mode='restore',versionname='before_calspectralflags')
Next we apply the calibration solutions to each source individually, using the gainfield parameter to specify which calibrators solutions should be applied from each of the gaintable calibration tables.
3C279 bandpass calibrator
# In CASA
applycal(vis='Band7multi_april22.ms',field='0',
gaintable=['bandpass.bcal','intphase.gcal','flux.cal'],
interp=['nearest','nearest','nearest'],
gainfield=['0','0','0'],flagbackup=T)
Titan absolute flux calibrator
# In CASA
applycal(vis='Band7multi_april22.ms',field='1',
gaintable=['bandpass.bcal','intphase.gcal','flux.cal'],
interp=['nearest','nearest','nearest'],
gainfield=['0','1','1'],flagbackup=T)
For the secondary phase calibrator, we apply the solutions from the primary phase calibrator. This will allow us to check how well the phase transfer to the science target (TW Hya) has worked.
# In CASA
applycal(vis='Band7multi_april22.ms',field='3',
gaintable=['bandpass.bcal','scanphase.gcal','flux.cal'],
interp=['nearest','nearest','nearest'],
gainfield=['0','4','4'],flagbackup=T)
Primary phase calibrator
# In CASA
applycal(vis='Band7multi_april22.ms',field='4',
gaintable=['bandpass.bcal','intphase.gcal','flux.cal'],
interp=['nearest','nearest','nearest'],
gainfield=['0','4','4'],flagbackup=T)
Science Target TW Hya. We apply the phase solutions from both the primary and secondary phase calibrators.
# In CASA
applycal(vis='Band7multi_april22.ms',field='2',
interp=['nearest','linear','linear'],
gaintable=['bandpass.bcal','scanphase.gcal','flux.cal'],
gainfield=['0','3,4','3,4'],flagbackup=T)
Now we can check the results of applying the calibration solutions.
# In CASA
plotms(vis='Band7multi_april22.ms',spw='',xaxis='time',yaxis='amp',field='',avgchannel='3840',
coloraxis='field',iteraxis='spw',ydatacolumn='corrected')
Looks good but there are some noisy points toward the end of the first scheduling block, mostly only affects the secondary phase calibrator (green). Lets check whether it improves when we apply this calibrators solutions to itself rather than transfering from the primary phase calibrator.
# In CASA
applycal(vis='Band7multi_april22.ms',field='3',
gaintable=['bandpass.bcal','intphase.gcal','flux.cal'],
interp=['nearest','nearest','nearest'],
gainfield=['0','3','3'],flagbackup=T)
# In CASA
plotms(vis='Band7multi_april22.ms',spw='',xaxis='time',yaxis='amp',field='',avgchannel='3840',
coloraxis='field',iteraxis='spw',ydatacolumn='corrected')
Indeed it looks better now. This is an indication that the amplitude calibration is not working well in this timerange. Since the effect is pretty small and we know we can self-cal later we let this go for now. Or if you want you can delve deeper into the data and try to solve the mystery...
Lets have a look at the corrected phase of the calibrators. If all has worked well, they should be near zero.
# In CASA
plotms(vis='Band7multi_april22.ms',spw='',xaxis='time',yaxis='phase',field='0,1,3,4',
avgchannel='3840',coloraxis='field',iteraxis='spw',ydatacolumn='corrected')
Indeed, these look very reasonable.
Lets check the spectral calibration now. First 3C279 and Titan
# In CASA
plotms(vis='Band7multi_april22.ms',spw='',xaxis='frequency',yaxis='amp',field='0,1',
avgtime='1e8',avgscan=T,coloraxis='field',ydatacolumn='corrected',
iteraxis='spw',xselfscale=T)
Flip through spws. These look good. Next the two phase calibrators.
# In CASA
plotms(vis='Band7multi_april22.ms',spw='',xaxis='frequency',yaxis='amp',field='3,4',
avgtime='1e8',avgscan=T,coloraxis='field',ydatacolumn='corrected',
iteraxis='spw',xselfscale=T)
Flip through spws. Also looks good, the secondary phase calibrator is a bit noisier but this is normal. The upswing at the lower frequency end of spw=3 will probably need to be flagged on the science target as well. It corresponds to the edge of an atmospheric feature.
Now for the exciting part, what does the science target look like...
# In CASA
plotms(vis='Band7multi_april22.ms',spw='',xaxis='frequency',yaxis='amp',field='2',
avgtime='1e8',avgscan=T,coloraxis='field',ydatacolumn='corrected',
iteraxis='spw',xselfscale=T)
We see HCO+(4-3) in spw=0 and CO(3-2) in spw=2 as expected.
Now check that the velocity of TWHydra is correct, it should be around 2.88 km/s
- CO 3-2 rest freq 345.79599 GHz
- HCO+ 4-3 rest freq 356.7342 GHz
# In CASA
plotms(vis='Band7multi_april22.ms',spw='2',xaxis='velocity',yaxis='amp',field='2',
avgtime='1e8',avgscan=T,transform=T,freqframe='LSRK',restfreq='345.79599GHz',
plotrange=[-10,15,0,0],coloraxis='spw')
# In CASA
plotms(vis='Band7multi_april22.ms',spw='0',xaxis='velocity',yaxis='amp',field='2',
avgtime='1e8',avgscan=T,transform=T,freqframe='LSRK',restfreq='356.7342GHz',
plotrange=[-10,15,0,0],coloraxis='spw')
Both lines show up at the expected velocity.
Split Calibrated Data
Now we split off the Science Target and interesting parts of the calibrators for imaging.
# In CASA
split(vis='Band7multi_april22.ms',outputvis='TWHydra_corrected.ms',
datacolumn='corrected',spw='',field='2')
# In CASA
split(vis='Band7multi_april22.ms',outputvis='J1037_corrected.ms',
datacolumn='corrected',spw='',field='4')
# In CASA
split(vis='Band7multi_april22.ms',outputvis='J1147_corrected.ms',
datacolumn='corrected',spw='',field='3')
So we check in more detail the effect of correcting the secondary phase calibrator with the primary, we'll redo the applycal and split a version with this calibration.
# In CASA
applycal(vis='Band7multi_april22.ms',field='3',
gaintable=['bandpass.bcal','scanphase.gcal','flux.cal'],
interp=['nearest','linear','linear'],
gainfield=['0','4','4'],flagbackup=T)
# In CASA
split(vis='Band7multi_april22.ms',outputvis='J1147_corrected_with4.ms',
datacolumn='corrected',spw='',field='3')
# In CASA
split(vis='Band7multi_april22.ms',outputvis='Titan_cont.ms',
datacolumn='corrected',spw='0',field='1')
# In CASA
split(vis='Band7multi_april22.ms',outputvis='Titan_CO3_2.ms',
datacolumn='corrected',spw='2',field='1')
# In CASA
split(vis='Band7multi_april22.ms',outputvis='Titan_unknown.ms',
datacolumn='corrected',spw='3',field='1')
# In CASA
split(vis='Band7multi_april22.ms',outputvis='3C279_CO3_2.ms',
datacolumn='corrected',spw='2:1900~2000',field='0')
Imaging Calibrators
This section will be added soon.
