IRAS16293 Band9 - Calibration for CASA 3.4: Difference between revisions

From CASA Guides
Jump to navigationJump to search
No edit summary
 
(40 intermediate revisions by 2 users not shown)
Line 5: Line 5:
*'''Details of the ALMA observations are provided at [[IRAS16293Band9]]
*'''Details of the ALMA observations are provided at [[IRAS16293Band9]]


*'''This portion of the guide covers calibration of the raw visibility data. To skip to the imaging portion of the guide, see: [[IRAS16293_Band9_-_Imaging]]'''.
*'''This portion of the guide covers calibration of the raw visibility data. To skip to the imaging portion of the guide, see: [[IRAS16293_Band9_-_Imaging_for_CASA_3.4]]'''.


<pre style="background-color: #ffa07a;">
WARNING: On June 15, 2012 the calibration guide and the final data products (calibrated science
data: IRAS16293_Band9_CalibratedMS_FIXED.tgz and reference images: IRAS16293_Band9_ReferenceImages_FIXED.tgz))
were changed to correct for a 1.2" position error in the phase calibrator (1625-254). Without
correction, the science images will suffer from a similar offset.
</pre>


==Overview==
==Overview==
Line 39: Line 45:
==Confirm your version of CASA==
==Confirm your version of CASA==


This guide has been written for CASA release 3.4.0.  Please confirm your version before proceeding. This guide will not work in earlier releases (<r19407).
This guide has been written for CASA release 3.4.0.  Please confirm your version before proceeding. This guide will not work in earlier versions (< r19874).
<source lang="python">
<source lang="python">
# In CASA
# In CASA
version = casalog.version()
version = casalog.version()
print "You are using " + version
print "You are using " + version
if (int(version.split()[4][1:-1]) < 19407):
if (int(version.split()[4][1:-1]) < 19874):
     print "\033[91m YOUR VERSION OF CASA IS TOO OLD FOR THIS GUIDE."
     print "\033[91m YOUR VERSION OF CASA IS TOO OLD FOR THIS GUIDE."
     print "\033[91m PLEASE UPDATE IT BEFORE PROCEEDING."
     print "\033[91m PLEASE UPDATE IT BEFORE PROCEEDING."
Line 57: Line 63:
http://casaguides.nrao.edu/index.php?title=Analysis_Utilities
http://casaguides.nrao.edu/index.php?title=Analysis_Utilities


for a full description and download instructions. If you do not wish to do this, see a CASA 3.3 version of one of the other ALMA guides for alternative (but slow) plotting options. Analysis Utilities are updated frequently so if its been a while since you installed it, its probably worth doing it again. If you are at an ALMA site or ARC, the analysis utilities are probably already installed and up to date.
for a full description and download instructions. If you do not wish to do this, see a CASA 3.3 version of one of the other ALMA guides for alternative (but slow) plotting options. Analysis Utilities are updated frequently so if its been a while since you installed it, it's probably worth doing it again. If you are at an ALMA site or ARC, the analysis utilities are probably already installed and up to date.
 


==Initial Inspection==
==Initial Inspection==


The first step we will do through all the calibration process is to define an array with the uid's that correspond to the datasets names. This will allow to make the calibration of the four datasets one after another, using a for-loop inside python. We will then calibrate the data individually and concatenate them at the end, before proceeding with the imaging part.
The first step we will do through all the calibration process is to define an array with the uid's that correspond to the dataset names. This will allow us to make the calibration of the four datasets one after another, using a for-loop inside python. We will then calibrate the data individually and concatenate them at the end, before proceeding with the imaging part.


Note that if you exit CASA and want to continue with the calibration using these arrays, you will have to re-issue the command again to make it available for the current CASA execution.
Note that if you exit CASA and want to continue with the calibration using these arrays, you will have to re-issue the command again to make it available for the current CASA execution.


To start, and give an example of this process, we will create txt format files for the output of the listobs task, which will give us useful information about the observations.
To start, and give an example of this process, we will create txt format files for the output of the {{listobs}} task, which will give us useful information about the observations.


<source lang="python">
<source lang="python">
Line 80: Line 85:
</source>
</source>


Note that after cutting and pasting a for-loop you often have to press return several times to make the command executes. The output will be sent to the CASA logger. Next there is an example of a useful part of the output that the first listobs of the previews command produces.
Note that after cutting and pasting a for-loop you often have to press return several times to make the command execute. The output will be sent to the CASA logger. Next, we show an example of a useful part of the output that the first listobs command produces.


<pre style="background-color: #fffacd;">
<pre style="background-color: #fffacd;">
Line 144: Line 149:
</pre>
</pre>


In the previous output you can see the ID that is assigned to each source, starting with the number 0. 1924-292 and 3c279 are the calibrators for bandpass, Juno for amplitud (flux), 1625-254 is our phase calibrator and nrao530 ph serves as a check (phase calibrator) source. The remaining 7 fields of IRAS16293-2422-a are the seven pointings for our mosaic of the target source.
In the previous output you can see the ID that is assigned to each source, starting with the number 0. 1924-292 and 3c279 are the calibrators for bandpass, Juno for amplitude (flux), 1625-254 is our phase calibrator and nrao530 ph serves as a check (phase calibrator) source. The remaining 7 fields of IRAS16293-2422-a are the seven pointings for our mosaic of the target source.


Spectral windows are also marked with numbers from 0 to 24, with number 0 containing WVR information. Spws 17, 19, 21, and 23 contain the sience data (FDM mode). The CO (6-5) line emission is contained in spw 19. Spw 18, 20, 22, and 24 contain channel averages of the data from spectral windows 17, 19, 21, 23, respectively. These spws will not be used for the offline data reduction. All the remaining spw that appear in the section of "Sources", and that do not appear in the "Spectral Windows" section are related to WVR measurements for each antenna, so you will not need them for the calibration neither. Spws 9, 11, 13, and 15 are associated with tsys measurements, and we will apply these information to the science spws later.
Spectral windows are also marked with numbers from 0 to 24, with number 0 containing WVR information. Spws 17, 19, 21, and 23 contain the science data (FDM mode). The CO (6-5) line emission is contained in spw 19. Spw 18, 20, 22, and 24 contain channel averages of the data from spectral windows 17, 19, 21, 23, respectively. These spws will not be used for the offline data reduction. All the remaining spws that appear in the section of "Sources", and that do not appear in the "Spectral Windows" section are related to WVR measurements for each antenna, so you will not need them for the calibration either. Spws 9, 11, 13, and 15 are associated with tsys measurements, and we will apply these measurements to the science spws later.


Finally before we go further we explicitly save the current flag state of the data. If you ever decide to start completely over, you should substitute 'restore' for 'save' in the command below to restore the flag state to its original value.
Finally before we go further we explicitly save the current flag state of the data. If you ever decide to start completely over, you should substitute 'restore' for 'save' in the command below to restore the flag state to its original value.
Line 156: Line 161:
   flagmanager(vis=vis,mode='save',versionname='OriginalFlagState')
   flagmanager(vis=vis,mode='save',versionname='OriginalFlagState')
</source>
</source>


== Generation and visualization of the antenna positions, Tsys and WVR tables ==
== Generation and visualization of the antenna positions, Tsys and WVR tables ==


Next we will generate the tables we need to apply to the data, antenna position, system temperature and water vapor radiometer. Once the tables are generated we will produce plots of them and inspect them to make sure whether they have issues that might affect their application to the data. Whenever we see an odd behavior in the tables we need to flag the corresponding science data to prevent wrong results in the calibration steps.  
Next we will generate the tables we need to apply to the data: antenna position, system temperature and water vapor radiometer. Once the tables are generated we will produce plots of them and inspect them to make sure whether they have issues that might affect their application to the data. Whenever we see an odd behavior in the tables we need to flag the corresponding science data to prevent wrong results in the calibration steps.  


=== System Temperature ===
=== System Temperature ===
Line 173: Line 177:
</source>
</source>


The next command, that comes from the Analysis Utils package will plot the tsys in the next way: it will produce many plots, each one of them will show an antenna, with the four spw that tsys tables cover, for all the targets, and with different colors for different times, so you can trace the behavior for tsys with time, among others.
The next command, that comes from the Analysis Utils package will plot the tsys in the next way: it will produce many plots, each one of them will show an antenna, with the four spws that the tsys tables cover, for all the targets, and with different colors for different times, so you can trace the behavior for tsys with time, among others.
Note that in spw 19, the overlap with the tsys spw (11) is not set correctly. This is due to an error in the frequencies for the tsys when the observations were done. You do not have to worry about this, since any issue coming from that error has already been fixed. Note, however that the portions of the spectra that do not have tsys information cannot be used. This does not represent a problem, since that part corresponds to the edge of the baseband. Also note that the CO (6-5) line is not affected by this.
Note that in spw 19, the overlap with the tsys spw (11) is not set correctly. This is due to an error in the frequencies for the tsys when the observations were done. You do not have to worry about this, since any issue coming from that error has already been fixed. Note, however that the portions of the spectra that do not have tsys information cannot be used. This does not represent a problem, since that part corresponds to the edge of the baseband. Also note that the CO (6-5) line is not affected by this.
In Figure 1 you will see the corresponding plot for one of the datasets (X90c) showing antenna 0 (DA41).
In Figure 1 you will see the corresponding plot for one of the datasets (X90c) showing antenna 0 (DA41).
Line 193: Line 197:
Another useful set of plots is generated by the next command. Each resulting plot shows all the antennas in a scan, for one spectral windows. This is useful to check Tsys for all the antennas in a certain time. In Figure 2 you can see an example.
Another useful set of plots is generated by the next command. Each resulting plot shows all the antennas in a scan, for one spectral windows. This is useful to check Tsys for all the antennas in a certain time. In Figure 2 you can see an example.


[[File:uid___A002_X3d55cb_X90c.ms.antenna.tdm.tsys.spw9.t1.png|200px|thumb|right|'''Fig. 2.''' Second sample plot for tsys plotting. The first four scans of the observations for spw 17 are showed, with all the antennas per plot.]]
[[File:uid___A002_X3d55cb_X90c.ms.antenna.tdm.tsys.spw9.t1.png|200px|thumb|right|'''Fig. 2.''' Second sample plot for tsys plotting. The first four scans of the observations for spw 17 are shown, with all the antennas per plot.]]


<source lang="python">
<source lang="python">
Line 223: Line 227:
  X39b
  X39b
  500 to 800 56 to 65%
  500 to 800 56 to 65%
  DA43 Tsys crazy for spw=23
  DA43 Tsys extremely high for spw=23
  DV05 ripples all spw and one time with bad YY  
  DV05 ripples all spw and one time with bad YY  
  DV18 crazy for spw=23
  DV18 extremely high for spw=23
</pre>
</pre>


Line 254: Line 258:
Note that in all datasets, DV15 has bad wvr behavior.
Note that in all datasets, DV15 has bad wvr behavior.


[[File:uid___A002_X3d55cb_X90c.ms.wvr.smooth.2.png|200px|thumb|right|'''Fig. 3.''' Phase corrections as funcion of time for the dataset X90c, where is shown an example of the odd behavior of DV15-related baselines.]]
[[File:uid___A002_X3d55cb_X90c.ms.wvr.smooth.2.png|200px|thumb|right|'''Fig. 3.''' Phase corrections as function of time for the dataset X90c, where you can see an example of the odd behavior of DV15-related baselines.]]


<source lang="python">
<source lang="python">
Line 274: Line 278:


# For the data taken on April 16
# For the data taken on April 16
antennas='DV18,DV12,DV09,DV10,DV13,DV05,DA41,DV14,DA43,DV17,DV07,DV15'
parameter=[0.00684754783288,-0.0188271608204,-0.00732900947332,
0.000481227821554,-2.48519197708e-05,-0.00048066949752,
0.000741523189806,-0.000843518709779,-0.000409113120338,
0.00111751207514,-0.000891124202747,-0.000492116263895,
0.00100313175518,-0.00084909900337,-0.000594194887412,
0.00105206234994,-0.000857719516042,-0.00059565765641,
0.00060776527971,-0.000525096431375,-0.000386487226933,
0.000521688542419,-0.000144920371688,-0.000258341151538,
0.00104740676978,-0.000899314162014,-0.000483014592547,
0.000920921096116,-0.00106288533016,-0.000546355296173,
0.000867350985729,-0.000990079062422,-0.00038123971161,
0.000567245762795,-0.000422531738877,-0.000741700641811]
gencal(vis = 'uid___A002_X3d4118_X39b.ms',
gencal(vis = 'uid___A002_X3d4118_X39b.ms',
   caltable = 'uid___A002_X3d4118_X39b.ms.antpos',
   caltable = 'uid___A002_X3d4118_X39b.ms.antpos',
   caltype = 'antpos',
   caltype = 'antpos',antenna=antennas,parameter=parameter)
  antenna = 'DV18,DV12,DV09,DV10,DV13,DV05,DA41,DV14,DA43,DV17,DV07,DV15',
  parameter = [0.00684754783288,-0.0188271608204,-0.00732900947332,0.000481227821554,-2.48519197708e-05,-0.00048066949752,0.000741523189806,-0.000843518709779,-0.000409113120338,0.00111751207514,-0.000891124202747,-0.000492116263895,0.00100313175518,-0.00084909900337,-0.000594194887412,0.00105206234994,-0.000857719516042,-0.00059565765641,0.00060776527971,-0.000525096431375,-0.000386487226933,0.000521688542419,-0.000144920371688,-0.000258341151538,0.00104740676978,-0.000899314162014,-0.000483014592547,0.000920921096116,-0.00106288533016,-0.000546355296173,0.000867350985729,-0.000990079062422,-0.00038123971161,0.000567245762795,-0.000422531738877,-0.000741700641811])


# and for April 17
# and for April 17
 
antennas='DV05,DV09,DV10,DV13,DV12,DA41,DV14,DA43,DV17,DV07,DV15'
parameter=[0.00105206234994,-0.000857719516042,-0.00059565765641,
0.000741523189806,-0.000843518709779,-0.000409113120338,
0.00111751207514,-0.000891124202747,-0.000492116263895,
0.00100313175518,-0.00084909900337,-0.000594194887412,
0.000481227821554,-2.48519197708e-05,-0.00048066949752,
0.00060776527971,-0.000525096431375,-0.000386487226933,
0.000521688542419,-0.000144920371688,-0.000258341151538,
0.00104740676978,-0.000899314162014,-0.000483014592547,
0.000920921096116,-0.00106288533016,-0.000546355296173,
0.000867350985729,-0.000990079062422,-0.00038123971161,
0.000567245762795,-0.000422531738877,-0.000741700641811]
for vis in range(1,4):
for vis in range(1,4):
   gencal(vis=rawdata[vis],
   gencal(vis=rawdata[vis],
     caltable=rawdata[vis]+'.antpos',
     caltable=rawdata[vis]+'.antpos',
     caltype = 'antpos',
     caltype ='antpos',antenna=antennas,parameter=parameter)
    antenna = 'DV05,DV09,DV10,DV13,DV12,DA41,DV14,DA43,DV17,DV07,DV15',
    parameter = [0.00105206234994,-0.000857719516042,-0.00059565765641,0.000741523189806,-0.000843518709779,-0.000409113120338,0.00111751207514,-0.000891124202747,-0.000492116263895,0.00100313175518,-0.00084909900337,-0.000594194887412,0.000481227821554,-2.48519197708e-05,-0.00048066949752,0.00060776527971,-0.000525096431375,-0.000386487226933,0.000521688542419,-0.000144920371688,-0.000258341151538,0.00104740676978,-0.000899314162014,-0.000483014592547,0.000920921096116,-0.00106288533016,-0.000546355296173,0.000867350985729,-0.000990079062422,-0.00038123971161,0.000567245762795,-0.000422531738877,-0.000741700641811])
</source>
</source>


Line 306: Line 330:
flagdata(vis=vis,antenna='DA43,DV18',spw='23',flagbackup=F)
flagdata(vis=vis,antenna='DA43,DV18',spw='23',flagbackup=F)


# For now leave DV05 in but keep an eye on it.
</source>
</source>


== Applying antpos, tsys, and wvr tables and splitting the data ==
== Applying antpos, tsys, and wvr tables and splitting the data ==


To apply the Tsys tables we need to separate sources with Tsys measurements of their own from those that
To apply the Tsys tables we need to separate sources that have Tsys measurements of their own from those that do not. Since only one IRAS16293 field (id=4) has Tsys, it goes into
do not. Since only one IRAS16293 field (id=4) has Tsys, it goes into
the "not" category. As you probably noted from the messages from
the "not" category. As you probably noted from the messages from
plotbandpass, the Tsys measurements on Juno did not yield usable
plotbandpass, the Tsys measurements on Juno did not yield usable
Line 336: Line 357:
</source>
</source>


To apply the Tsys we need to explicitly tell applycal which TDM (128
To apply the Tsys we need to explicitly tell {{applycal}} which TDM (128
channels) spws go with which FDM (3840 channels). An entry is needed
channels) spws go with which FDM spw (3840 channels). An entry is needed
for every spw, because the position in the list corresponds to the
for every spw, because the position in the list corresponds to the
spw id to be corrected.
spw id to be corrected.
Line 353: Line 374:
The helper gets a little confused in this case with spws 19 and 20
The helper gets a little confused in this case with spws 19 and 20
because the FDM was not perfectly positioned within the TDM window
because the FDM was not perfectly positioned within the TDM window
as you can see from the Tsys plots. In a case like this we need to
as you can see from the Tsys plots. In a case like this, we need to
put in by hand the TDM window that should match the FDM window.
put in by hand the TDM window that should match the FDM window.


As you could see from the first run of listobs, 1924-292 is the bandpass calibrator for three of the datasets, and 3c279 is for one of them. For this reason, the application of the tables is split in two parts.
As you could see from the first run of listobs, 1924-292 is the bandpass calibrator for three of the datasets, and 3c279 is for one of them. For this reason, the application of the tables is split into two parts.


For the {{aplycal}} task, note that the first value in the interp parameter for each gaintable gives the
For the {{applycal}} task, note that the first value in the interp parameter for each gaintable gives the desired time interpolation type.  The second parameter indicates the
desired time interpolation type.  The second parameter indicates the
desired frequency axis interpolation. If the second value is not set
desired frequency axis interpolation. If the second value is not set
it is assumed to be linear, but only if the input table has a
it is assumed to be linear, but only if the input table has a
Line 430: Line 450:
</source>
</source>


== Fix Phase Calibrator Position ==
The position used for the phase calibrator (1625-254) in these observations is offset by about 1.2"
toward postive R.A from its correct position. Below we correct the data for this offset. It is best to do this step before doing any of the gain calibration. It is notable that simply shifting the images for this offset will not give as correct a result due to phase delays caused by sky rotation during the course of each observation.
<source lang="python">
data=['uid___A002_X3d4118_X39b.antwvrtsys.ms',
      'uid___A002_X3d55cb_X575.antwvrtsys.ms',
      'uid___A002_X3d55cb_X90c.antwvrtsys.ms',
      'uid___A002_X3d55cb_Xb50.antwvrtsys.ms']
</source>
List the data to see current position of 1625-254
<source lang="python">
for vis in data:
  listobs(vis=vis,listfile=vis+'.listobs',verbose=True)
</source>
Incorrect position from listobs: 1625-254  16h25m46.98000s -25d27m38.3300s J2000
From the [http://www.vla.nrao.edu/astro/calib/manual/csource.html
EVLA calibrator manual] we find that the correct position is: 16h25m46.891639s -25d27m38.326880s J2000
The task fixvis still requires scratch columns, so we must create
these first (this takes a while).
<source lang="python">
for vis in data:
  clearcal(vis=vis)
</source>
Next use the task fixvis to correct the position of the calibrator in
the data and the header. It will also recalculate UVWs, but this
correction is very small for a 1.2" shift.
<source lang="python">
for vis in data:
  fixvis(vis=vis,outputvis=vis+'.fixed',field='1625-254',reuse=False,
        phasecenter='J2000 16h25m46.891639s -25d27m38.326880s')
</source>
Define the new dataset names that we will use here on out in the calibration.
<source lang="python">
data=['uid___A002_X3d4118_X39b.antwvrtsys.ms.fixed',
      'uid___A002_X3d55cb_X575.antwvrtsys.ms.fixed',
      'uid___A002_X3d55cb_X90c.antwvrtsys.ms.fixed',
      'uid___A002_X3d55cb_Xb50.antwvrtsys.ms.fixed']
</source>
Next, we need to remove the model data column, as fixvis currently
rotates both the data AND the model data column. If you do not do the
following step, you will double correct for the position offset.
<source lang="python">
for vis in data:
  delmod(vis=vis,scr=True,otf=False)
</source>
Check the new position in listobs
<source lang="python">
for vis in data:
  listobs(vis=vis,listfile=vis+'listobs_afterfix',verbose=F)
</source>


== Data inspection ==
== Data inspection ==


We now need to check any bad behavior in the data through several plots. Once problems are identified, data can be flagged. But before that, we need to run again {{listobs}} to check that the split worked as expected. We will define our new array of split datasets, along with a list of intents that will be useful in the next steps.
We now need to check for any bad behavior in the data through several plots. Once problems are identified, data can be flagged. But before that, we need to run again {{listobs}} to check that the split worked as expected. We will define our new array of split datasets, along with a list of intents that will be useful in the next steps.




<source lang="python">
<source lang="python">
# New array of datasets  
# New array of datasets  
data=['uid___A002_X3d4118_X39b.antwvrtsys.ms',
data=['uid___A002_X3d4118_X39b.antwvrtsys.ms.fixed',
       'uid___A002_X3d55cb_X575.antwvrtsys.ms',
       'uid___A002_X3d55cb_X575.antwvrtsys.ms.fixed',
       'uid___A002_X3d55cb_X90c.antwvrtsys.ms',
       'uid___A002_X3d55cb_X90c.antwvrtsys.ms.fixed',
       'uid___A002_X3d55cb_Xb50.antwvrtsys.ms']
       'uid___A002_X3d55cb_Xb50.antwvrtsys.ms.fixed']


# Match up intents with source names, for the lists below, it is important  
# Match up intents with source names, for the lists below, it is important  
Line 523: Line 609:
In Figure 5 you can see the output of the following plotms command, which runs with X39b. By clicking the "Next" arrow in plotms you can access the remaining spws, since the command was executed with the option iteraxis='spw'.
In Figure 5 you can see the output of the following plotms command, which runs with X39b. By clicking the "Next" arrow in plotms you can access the remaining spws, since the command was executed with the option iteraxis='spw'.


[[File:uid___A002_X3d4118_X39b.antwvrtsys.ms.freqamp_spw0.png|200px|thumb|right|'''Fig. 5.''' {{plotms}} result for amplitude vs time for all the sources, which are displayed with different colors. The plot shows spw 0 for the dataset X39b.]]
[[File:uid___A002_X3d4118_X39b.antwvrtsys.ms.fixed.time.spw0.png|200px|thumb|right|'''Fig. 5.''' {{plotms}} result for amplitude vs time for all the sources, which are displayed with different colors. The plot shows spw 0 for the dataset X39b.]]


<source lang="python">
<source lang="python">
Line 533: Line 619:
       iteraxis='spw',ydatacolumn='data',yselfscale=True)
       iteraxis='spw',ydatacolumn='data',yselfscale=True)
</source>
</source>
For the next plotms commands, inspect each dataset, noting any problems that you notice.
For the next set of plotms commands, inspect each dataset, noting any problems that you notice.


<source lang="python">
<source lang="python">
Line 559: Line 645:
       iteraxis='baseline',ydatacolumn='data',yselfscale=True)
       iteraxis='baseline',ydatacolumn='data',yselfscale=True)
</source>
</source>


== Flagging ==
== Flagging ==
Line 566: Line 651:


<source lang="python">
<source lang="python">
data=['uid___A002_X3d4118_X39b.antwvrtsys.ms',
data=['uid___A002_X3d4118_X39b.antwvrtsys.ms.fixed',
       'uid___A002_X3d55cb_X575.antwvrtsys.ms',
       'uid___A002_X3d55cb_X575.antwvrtsys.ms.fixed',
       'uid___A002_X3d55cb_X90c.antwvrtsys.ms',
       'uid___A002_X3d55cb_X90c.antwvrtsys.ms.fixed',
       'uid___A002_X3d55cb_Xb50.antwvrtsys.ms']
       'uid___A002_X3d55cb_Xb50.antwvrtsys.ms.fixed']


# Back up flag state in case you want to start over.
# Back up flag state in case you want to start over.
Line 580: Line 665:
</source>
</source>


Here is the flag that we will do. If you note some other data that needs to be flagged, proceed with that.
Here is the flagging that we will do. If you note some other data that needs to be flagged, proceed with that as well.


<source lang="python">
<source lang="python">
# DV07 low amp after certain time. DV17 also but not as bad.
# DV07 low amp after certain time. DV17 also but not as bad.
flagdata(vis='uid___A002_X3d55cb_X575.antwvrtsys.ms',
flagdata(vis='uid___A002_X3d55cb_X575.antwvrtsys.ms.fixed',
         antenna='DV07', timerange='>05:20:00',flagbackup=F)
         antenna='DV07', timerange='>05:20:00',flagbackup=F)
</source>
</source>


Now, during the procedure of calibration some problems in the data will show up. Those data needed to be flagged and do the calibration again. In order to save you time with this, we note here those problems, so you can flag the data now.
Now, during the procedure of calibration some problems in the data will show up. Those data need to be flagged and then the calibration should be repeated. In order to save you time with this, we note here those problems, so you can flag the data now.


<source lang="python">
<source lang="python">
# AFTER INITIAL CALIBRATION INSPECTION
# PROBLEMS DISCOVERED AFTER INITIAL CALIBRATION INSPECTION


# flag low elevation scans on 1625-254 and IRAS16293
# flag low elevation scans on 1625-254 and IRAS16293
flagdata(vis='uid___A002_X3d55cb_Xb50.antwvrtsys.ms',
flagdata(vis='uid___A002_X3d55cb_Xb50.antwvrtsys.ms.fixed',
         timerange='>11:54:00', field='',flagbackup=F)
         timerange='>11:54:00', field='',flagbackup=F)


# flag low gains on DV02 on Juno
# flag low gains on DV02 on Juno
flagdata(vis='uid___A002_X3d55cb_Xb50.antwvrtsys.ms',
flagdata(vis='uid___A002_X3d55cb_Xb50.antwvrtsys.ms.fixed',
         antenna='DV02', field='Juno',flagbackup=F)
         antenna='DV02', field='Juno',flagbackup=F)


flagdata(vis='uid___A002_X3d4118_X39b.antwvrtsys.ms',
flagdata(vis='uid___A002_X3d4118_X39b.antwvrtsys.ms.fixed',
         antenna='DV02', field='Juno',flagbackup=F)
         antenna='DV02', field='Juno',flagbackup=F)
</source>
</source>


== Calibration ==
== Calibration ==


Now we can start with the calibration itself. First, we will perform the bandpass calibration using 1924-292 and 3c279. As before, we define our list of data and match the sources with intents. Also, we will set our reference antenna (one close to the center of the array and without problems, like delays). As you can see, two different intervals for channels are used, which make sense in a bit.
Now we can start with the calibration itself. First, we will perform the bandpass calibration using 1924-292 and 3c279. As before, we define our list of data and match the sources with intents. Also, we will set our reference antenna (one close to the center of the array and without problems, like delays). As you can see, two different intervals for channels are used, which will make sense in a bit.


<source lang="python">
<source lang="python">
data=['uid___A002_X3d4118_X39b.antwvrtsys.ms',
data=['uid___A002_X3d4118_X39b.antwvrtsys.ms.fixed',
       'uid___A002_X3d55cb_X575.antwvrtsys.ms',
       'uid___A002_X3d55cb_X575.antwvrtsys.ms.fixed',
       'uid___A002_X3d55cb_X90c.antwvrtsys.ms',
       'uid___A002_X3d55cb_X90c.antwvrtsys.ms.fixed',
       'uid___A002_X3d55cb_Xb50.antwvrtsys.ms']
       'uid___A002_X3d55cb_Xb50.antwvrtsys.ms.fixed']
# Match up intents with source names
# Match up intents with source names
pcal='1625-254'
pcal='1625-254'
Line 637: Line 721:
<source lang="python">
<source lang="python">
for vis in range(len(data)):
for vis in range(len(data)):
   delmod(vis=data[vis]),
   delmod(vis=data[vis])
   setjy(vis=data[vis],field=fluxcal,standard='Butler-JPL-Horizons 2010',
   setjy(vis=data[vis],field=fluxcal,standard='Butler-JPL-Horizons 2010',
         scalebychan=True,usescratch=False)
         scalebychan=True,usescratch=False)
</source>
</source>
Now we proceed with the bandpass, which will be done in two steps. In the first one, we will use the block of channels in the center of the spws, 1200~1500, to calculate gains in phase, using a solution interval of 30 seconds. This will give us the variation of the phase through all the observation. We will use this table later to execute the {{bandpass}} itself. Note that we are using a minimum of signal to noise ratio of 2 to accept the solutions. This is low but needed since the calibrators are not very strong in band 9.
Now we proceed with the bandpass calibration, which will be done in two steps. In the first one, we will use a small block of channels in the center of the spws, 1200~1500, to calculate gains in phase, using a solution interval of 30 seconds. This will give us the variation of the phase throughout the observation. We will use this table later to execute the {{bandpass}} itself. Note that we are using a minimum signal to noise ratio of 2 to accept the solutions. This is low but necessary since the calibrators are not very strong in band 9.


<source lang="python">
<source lang="python">
Line 651: Line 735:
</source>
</source>


Next, we set our bandpass command to use our previously generated gain tables. Since we do not have a high signal to noise ratio per channel, we use a polynomial option to calculate the solutions for the bandpass. In the {{bandpass}} task, the options degamp and degphase will set the maximum degree of the polynomial that the task can use to calculate solutions. In the work log you can note what is the actual degree that the task is using. Note that the combination of solint='inf' and combine='scan' will result in a solution per scan for the calibrator.
Next, we set our bandpass command to use our previously generated gain tables. Since we do not have a high signal to noise ratio per channel, we use a polynomial option to calculate the solutions for the bandpass. In the {{bandpass}} task, the options degamp and degphase will set the maximum degree of the polynomial that the task can use to calculate solutions. In the work log you can note what is the actual degree that the task is using. Note that the combination of solint='inf' and combine='scan' will result in one solution per scan for the calibrator.
 
[[File:uid___A002_X3d55cb_X90c.antwvrtsys.ms.fixed.bandpass_bpoly.DV09.spw0.t1_34.png|200px|thumb|right|'''Fig. 6.''' Bandpass solutions for the dataset X90c (all 4 spws) showing antenna DV09. Both amplitude and phases solutions are plotted. These data use the weaker bandpass calibrator J1924-224.]]
 
[[File:uid___A002_X3d55cb_X575.antwvrtsys.ms.fixed.bandpass_bpoly.DV09.spw0.t1.png|200px|thumb|right|'''Fig. 7.''' Bandpass solutions for the dataset X575 (all 4 spws) showing antenna DV09. Both amplitude and phases solutions are plotted. These data use the stronger bandpass calibrator 3C279.]]
 


<source lang="python">
<source lang="python">
Line 670: Line 759:


To plot all our tables we will use both our AnalysisUtils and {{plotcal}}. The aU.plotbandpass command will create plot files for each combination of dataset and antenna, for both amplitude and phase. Inspect all these plots to make sure that the bandpass solutions look good.
To plot all our tables we will use both our AnalysisUtils and {{plotcal}}. The aU.plotbandpass command will create plot files for each combination of dataset and antenna, for both amplitude and phase. Inspect all these plots to make sure that the bandpass solutions look good.
In Figure 6 we show an output sample for this command.
In Figures 6 and 7 we show an output sample for this command for the two bandpass calibrators used -- note how noisy the weaker bandpass calibrator (J1924) is without BPOLY.
 
[[File:uid___A002_X3d55cb_X90c.antwvrtsys.ms.bandpass_bpoly.DA41.spw0.t1_34.png|200px|thumb|right|'''Fig. 6.''' Bandpass solutions for the dataset X90c (spw 0) showing antenna DA41. Both amplitude and phases solutions are plotted.]]


<source lang="python">
<source lang="python">
Line 690: Line 777:
</source>
</source>


The next step in the calibration is to calculate amplitude and phase gains vs time for our calibrators. The ideal case here is to have a solution per integration of the data, but in this case we will need to use a solution interval of 30 seconds to avoid having many failed solutions, especially in the week phase calibrator. First, we will calculate gains for phase and later using that information, we will solve for amplitude and phase. In all the next three executions we are using the bandpass calibration table, as it provides the gains for phase and amplitude vs frequency.
The next step in the calibration is to calculate amplitude and phase gains vs time for our calibrators. The ideal case here is to have a solution per integration of the data, but in this case we will need to use a solution interval of 30 seconds to avoid having many failed solutions, especially in the weak phase calibrator. First, we will calculate gains for phase and later using that information, we will solve for amplitude and phase. In all the next three executions we are using the bandpass calibration table, as it provides the gains for phase and amplitude vs frequency.


<source lang="python">
<source lang="python">
#Using 30s (5 integrations) per solution to avoid many failed solution of the week calibrator.
#Using 30s (5 integrations) per solution to avoid many failed solution of the weak calibrator.
for vis in range(len(data)):
for vis in range(len(data)):
   gaincal(vis=data[vis],caltable=data[vis]+'.intphase.gcal',
   gaincal(vis=data[vis],caltable=data[vis]+'.intphase.gcal',
Line 700: Line 787:
         gaintable=[data[vis]+'.bandpass_bpoly.bcal'])
         gaintable=[data[vis]+'.bandpass_bpoly.bcal'])
</source>
</source>
To plot our *.intphase.gcal tables we use the next command, which will create files containing the phase gains vs time for all the atennas and for all the datasets.
To plot our *.intphase.gcal tables we use the next command, which will create files containing the phase gains vs time for all the antennas and for all the datasets.
In Figure 7 we show an example of such plots. Again, you will need to check all plots to make sure the solutions are good.
In Figure 7 we show an example of such plots. Again, you will need to check all plots to make sure the solutions are good.


[[File:uid___A002_X3d55cb_X90c.antwvrtsys.ms.intphase.gcal.spw0.ant5_9.png|200px|thumb|right|'''Fig. 7.''' Phase gain solutions vs time for 5 antennas for all calibrators in the case of X90c.]]
[[File:uid___A002_X3d55cb_X90c.antwvrtsys.ms.fixed.intphase.gcal.spw0.ant5_9.png|200px|thumb|right|'''Fig. 8.''' Phase gain solutions vs time for 5 antennas for all calibrators in the case of X90c.]]


<source lang="python">
<source lang="python">
Line 730: Line 817:
The next command will produce many plots, and like the previous one, you will get one for each dataset and for each antenna. See example in Figure 8.
The next command will produce many plots, and like the previous one, you will get one for each dataset and for each antenna. See example in Figure 8.


[[File:uid___A002_X3d55cb_X90c.antwvrtsys.ms.scanphase.gcal.spw0.ant5_9.png|200px|thumb|right|'''Fig. 8.''' Phase gain solutions vs time for all the calibrators in the case of X90c showing the same 5 antennas from previous figure. Note that there is only one solution per scan for each source.]]
[[File:uid___A002_X3d55cb_X90c.antwvrtsys.ms.fixed.scanphase.gcal.spw0.ant5_9.png|200px|thumb|right|'''Fig. 9.''' Phase gain solutions vs time for all the calibrators in the case of X90c showing the same 5 antennas from previous figure. Note that there is only one solution per scan for each source.]]


<source lang="python">
<source lang="python">
Line 746: Line 833:
</source>
</source>


Finally, in the next {{gaincal}} we will solve for amplitude and phase. We ill use the gain phase calibration table produced before. We will get one solution per scan for all our calibrators.
Finally, in the next {{gaincal}} we will solve for amplitude and phase. We will use the gain phase calibration table produced before. We will get one solution per scan for all our calibrators.


<source lang="python">
<source lang="python">
Line 759: Line 846:
We now can check the resulting plots with the next plotcal executions. In Figure 9 we show an example of these plots.
We now can check the resulting plots with the next plotcal executions. In Figure 9 we show an example of these plots.


[[File:uid___A002_X3d55cb_X90c.antwvrtsys.ms.amp.gcal.spw0.ant5_9.png|200px|thumb|right|'''Fig. 9.''' Gain amplitudes for all the calibrators, as function of time. A plot for X90c is shown with 4 antennas in it. As in the previous plot, there is only one solution per scan for each source.]]
[[File:uid___A002_X3d55cb_X90c.antwvrtsys.ms.fixed.amp.gcal.spw0.ant5_9.png|200px|thumb|right|'''Fig. 10.''' Gain amplitudes for all the calibrators, as function of time. A plot for X90c is shown with 5 antennas in it. As in the previous plot, there is only one solution per scan for each source.]]


<source lang="python">
<source lang="python">
Line 776: Line 863:
We can plot the same tables but in a different way that will allow us to look for higher abnormal gains in the solutions. With the next command we will get only four plots for each spw. Each of those plots shows all the gains for all the antennas for all the sources (See Figure 10 for an example). Below in the box, there are some comments to focus you on some data, so you can double check. Before continuing, make sure you check all the calibration tables.
We can plot the same tables but in a different way that will allow us to look for higher abnormal gains in the solutions. With the next command we will get only four plots for each spw. Each of those plots shows all the gains for all the antennas for all the sources (See Figure 10 for an example). Below in the box, there are some comments to focus you on some data, so you can double check. Before continuing, make sure you check all the calibration tables.


[[File:uid___A002_X3d55cb_X90c.antwvrtsys.ms.amp.png|200px|thumb|right|'''Fig. 10.''' Gain amplitudes for all the calibrators, as function of time (X90c) for the spw 2. The plot shows all the antennas in each plot. Check for higher gains that deviate from the average.]]
[[File:uid___A002_X3d55cb_X90c.antwvrtsys.ms.fixed.amp.png|200px|thumb|right|'''Fig. 11.''' Gain amplitudes for all the calibrators, as function of time (X90c), all antennas for each spw. Check for higher gains that deviate from the average.]]


<source lang="python">
<source lang="python">
Line 784: Line 871:
   plotcal(caltable=vis+'.amp.gcal',
   plotcal(caltable=vis+'.amp.gcal',
         xaxis='time',yaxis='amp',antenna='',field='',
         xaxis='time',yaxis='amp',antenna='',field='',
         iteration='spw',subplot=221,poln='',spw='',
         iteration='spw',subplot=411,poln='',spw='',
         showgui=False,figfile=vis+'.amp.png')
         showgui=False,figfile=vis+'.amp.png')
# X90c good
# X90c good
Line 792: Line 879:
</source>
</source>


Now that the gain calibration is done, we need to set the flux for our calibrators. For this we will use  
==== Set Absolute Flux Scale ====
Juno, our primary flux calibrator. We will do this by using {{fluxscale}}.  
 
Now that the gain calibration is done, we need to set the flux for our calibrators. For this we will use Juno, our primary flux calibrator. We will do this by using {{fluxscale}}.  
We then will transfer the flux information from our phase calibrator to our science target.
We then will transfer the flux information from our phase calibrator to our science target.
For these datasets, we note the following for this step:
For these datasets, we note the following for this step:
Line 808: Line 896:
</source>
</source>


Below we copy the fluxscale numbers from the logger window:
<pre style="background-color: #fffacd;">
#:calibrater::open    Opening MS: uid___A002_X3d4118_X39b.antwvrtsys.ms.fixed for calibration.
#:Calibrater:: Initializing nominal selection to the whole MS.
#:calibrater::fluxscale        Beginning fluxscale--(MSSelection version)-------
#:::    Assuming all non-reference fields are transfer fields.
#:::    Found reference field(s): Juno
#:::    Found transfer field(s):  1924-292 nrao530 ph 1625-254
#:::    Flux density for 1924-292 in SpW=0 is: 2.29951 +/- 0.0615168 (SNR = 37.3802, N = 24)
#:::    Flux density for 1924-292 in SpW=1 is: 2.3884 +/- 0.0706186 (SNR = 33.8212, N = 24)
#:::    Flux density for 1924-292 in SpW=2 is: 2.30261 +/- 0.0561279 (SNR = 41.0244, N = 24)
#:::    Flux density for 1924-292 in SpW=3 is: 2.27951 +/- 0.0682752 (SNR = 33.3871, N = 20)
#:::    Flux density for nrao530 ph in SpW=0 is: 0.627546 +/- 0.0366209 (SNR = 17.1363, N = 24)
#:::    Flux density for nrao530 ph in SpW=1 is: 0.64101 +/- 0.0406892 (SNR = 15.7538, N = 24)
#:::    Flux density for nrao530 ph in SpW=2 is: 0.615775 +/- 0.0324313 (SNR = 18.9871, N = 24)
#:::    Flux density for nrao530 ph in SpW=3 is: 0.607186 +/- 0.0407003 (SNR = 14.9184, N = 20)
#:::    Flux density for 1625-254 in SpW=0 is: 0.405301 +/- 0.0270957 (SNR = 14.9581, N = 24)
#:::    Flux density for 1625-254 in SpW=1 is: 0.414875 +/- 0.0274183 (SNR = 15.1313, N = 24)
#:::    Flux density for 1625-254 in SpW=2 is: 0.411922 +/- 0.0212975 (SNR = 19.3413, N = 24)
#:::    Flux density for 1625-254 in SpW=3 is: 0.406868 +/- 0.0302855 (SNR = 13.4344, N = 20)
#:Calibrater::fluxscale        Storing result in uid___A002_X3d4118_X39b.antwvrtsys.ms.fixed.flux.cal
#:::  Writing solutions to table: uid___A002_X3d4118_X39b.antwvrtsys.ms.fixed.flux.cal
#:calibrater::open    Opening MS: uid___A002_X3d55cb_X575.antwvrtsys.ms.fixed for calibration.
#:Calibrater:: Initializing nominal selection to the whole MS.
#:calibrater::fluxscale        Beginning fluxscale--(MSSelection version)-------
#:::    Assuming all non-reference fields are transfer fields.
#:::    Found reference field(s): Juno
#:::    Found transfer field(s):  3c279 nrao530 ph 1625-254
#:::    Flux density for 3c279 in SpW=0 is: 10.0511 +/- 0.170046 (SNR = 59.1085, N = 22)
#:::    Flux density for 3c279 in SpW=1 is: 11.0567 +/- 0.204718 (SNR = 54.0092, N = 22)
#:::    Flux density for 3c279 in SpW=2 is: 10.8403 +/- 0.190514 (SNR = 56.9002, N = 22)
#:::    Flux density for 3c279 in SpW=3 is: 10.9256 +/- 0.213188 (SNR = 51.2484, N = 22)
#:::    Flux density for nrao530 ph in SpW=0 is: 0.895265 +/- 0.0469496 (SNR = 19.0686, N = 22)
#:::    Flux density for nrao530 ph in SpW=1 is: 0.962985 +/- 0.057688 (SNR = 16.693, N = 22)
#:::    Flux density for nrao530 ph in SpW=2 is: 0.95216 +/- 0.063454 (SNR = 15.0055, N = 22)
#:::    Flux density for nrao530 ph in SpW=3 is: 0.953193 +/- 0.0667595 (SNR = 14.278, N = 22)
#:::    Flux density for 1625-254 in SpW=0 is: 0.633491 +/- 0.054662 (SNR = 11.5892, N = 22)
#:::    Flux density for 1625-254 in SpW=1 is: 0.635812 +/- 0.0564569 (SNR = 11.2619, N = 22)
#:::    Flux density for 1625-254 in SpW=2 is: 0.62294 +/- 0.0523949 (SNR = 11.8893, N = 22)
#:::    Flux density for 1625-254 in SpW=3 is: 0.640096 +/- 0.0646176 (SNR = 9.90592, N = 22)
#:Calibrater::fluxscale        Storing result in uid___A002_X3d55cb_X575.antwvrtsys.ms.fixed.flux.cal
#:::  Writing solutions to table: uid___A002_X3d55cb_X575.antwvrtsys.ms.fixed.flux.cal


Analyze all the numbers coming out in the log window and write them now, as they can serve you as a reference later on.
#:calibrater::open    Opening MS: uid___A002_X3d55cb_X90c.antwvrtsys.ms.fixed for calibration.
#:Calibrater:: Initializing nominal selection to the whole MS.
#:calibrater::fluxscale        Beginning fluxscale--(MSSelection version)-------
#:::    Assuming all non-reference fields are transfer fields.
#:::    Found reference field(s): Juno
#:::    Found transfer field(s):  1924-292 nrao530 ph 1625-254
#:::    Flux density for 1924-292 in SpW=0 is: 2.26482 +/- 0.0804822 (SNR = 28.1406, N = 22)
#:::    Flux density for 1924-292 in SpW=1 is: 2.30424 +/- 0.0688451 (SNR = 33.4698, N = 22)
#:::    Flux density for 1924-292 in SpW=2 is: 2.31636 +/- 0.063015 (SNR = 36.7589, N = 22)
#:::    Flux density for 1924-292 in SpW=3 is: 2.30091 +/- 0.0628643 (SNR = 36.6011, N = 22)
#:::    Flux density for nrao530 ph in SpW=0 is: 0.655383 +/- 0.0392507 (SNR = 16.6974, N = 22)
#:::    Flux density for nrao530 ph in SpW=1 is: 0.663148 +/- 0.0358838 (SNR = 18.4804, N = 22)
#:::    Flux density for nrao530 ph in SpW=2 is: 0.660628 +/- 0.0318434 (SNR = 20.7462, N = 22)
#:::    Flux density for nrao530 ph in SpW=3 is: 0.656224 +/- 0.0347546 (SNR = 18.8817, N = 22)
#:::    Flux density for 1625-254 in SpW=0 is: 0.423542 +/- 0.0340076 (SNR = 12.4543, N = 22)
#:::    Flux density for 1625-254 in SpW=1 is: 0.413005 +/- 0.0278793 (SNR = 14.8141, N = 22)
#:::    Flux density for 1625-254 in SpW=2 is: 0.438477 +/- 0.030868 (SNR = 14.2049, N = 22)
#:::    Flux density for 1625-254 in SpW=3 is: 0.434589 +/- 0.028559 (SNR = 15.2173, N = 22)
#:Calibrater::fluxscale        Storing result in uid___A002_X3d55cb_X90c.antwvrtsys.ms.fixed.flux.cal
#:::  Writing solutions to table: uid___A002_X3d55cb_X90c.antwvrtsys.ms.fixed.flux.cal


#:calibrater::open    Opening MS: uid___A002_X3d55cb_Xb50.antwvrtsys.ms.fixed for calibration.
#:Calibrater:: Initializing nominal selection to the whole MS.
#:calibrater::fluxscale        Beginning fluxscale--(MSSelection version)-------
#:::    Assuming all non-reference fields are transfer fields.
#:::    Found reference field(s): Juno
#:::    Found transfer field(s):  1924-292 nrao530 ph 1625-254
#:::    Flux density for 1924-292 in SpW=0 is: 2.39168 +/- 0.0655984 (SNR = 36.4594, N = 22)
#:::    Flux density for 1924-292 in SpW=1 is: 2.33502 +/- 0.0478833 (SNR = 48.7649, N = 22)
#:::    Flux density for 1924-292 in SpW=2 is: 2.26212 +/- 0.053858 (SNR = 42.0016, N = 22)
#:::    Flux density for 1924-292 in SpW=3 is: 2.28193 +/- 0.0551113 (SNR = 41.4059, N = 22)
#:::    Flux density for nrao530 ph in SpW=0 is: 0.686101 +/- 0.0515145 (SNR = 13.3186, N = 22)
#:::    Flux density for nrao530 ph in SpW=1 is: 0.696972 +/- 0.0392324 (SNR = 17.7652, N = 22)
#:::    Flux density for nrao530 ph in SpW=2 is: 0.683672 +/- 0.0422611 (SNR = 16.1773, N = 22)
#:::    Flux density for nrao530 ph in SpW=3 is: 0.69109 +/- 0.0405247 (SNR = 17.0535, N = 22)
#:::    Flux density for 1625-254 in SpW=0 is: 0.504527 +/- 0.0496797 (SNR = 10.1556, N = 22)
#:::    Flux density for 1625-254 in SpW=1 is: 0.482058 +/- 0.0293865 (SNR = 16.4041, N = 22)
#:::    Flux density for 1625-254 in SpW=2 is: 0.445422 +/- 0.0365173 (SNR = 12.1976, N = 22)
#:::    Flux density for 1625-254 in SpW=3 is: 0.458493 +/- 0.0345567 (SNR = 13.2678, N = 22)
#:Calibrater::fluxscale        Storing result in uid___A002_X3d55cb_Xb50.antwvrtsys.ms.fixed.flux.cal
#:::  Writing solutions to table: uid___A002_X3d55cb_Xb50.antwvrtsys.ms.fixed.flux.cal
</pre>
                                                                                                           
When derived fluxes are too high and nothing else appears wrong with
When derived fluxes are too high and nothing else appears wrong with
the data, the cause especially at Band 9 is likely decorrelation. So
the data, the cause (especially at Band 9) is likely to be decorrelation. So
we will favor the average of the lower values to explicitly set the
we will favor the average of the lower values to explicitly set the
flux densities based on the {{fluxscale}} results. We set the two
flux densities based on the {{fluxscale}} results. We set the flux of the two
bandpass calibrators for convenience of having fully calibrated
bandpass calibrators for the convenience of having a fully calibrated
dataset. The only ones that really matter are 1625 as the gain
dataset. The only ones that really matter are 1625 as the gain
calibrator and the check source nrao530.
calibrator and the check source nrao530.
Based on this, we proceed as following to make the changes that are needed:
Based on this, we proceed as follows to make the changes that are needed:


<source lang="python">
<source lang="python">
Line 827: Line 1,000:
flux3c279=[10.5,0,0,0]
flux3c279=[10.5,0,0,0]


datawith1924=['uid___A002_X3d4118_X39b.antwvrtsys.ms',
datawith1924=['uid___A002_X3d4118_X39b.antwvrtsys.ms.fixed',
       'uid___A002_X3d55cb_X90c.antwvrtsys.ms',
       'uid___A002_X3d55cb_X90c.antwvrtsys.ms.fixed',
       'uid___A002_X3d55cb_Xb50.antwvrtsys.ms']
       'uid___A002_X3d55cb_Xb50.antwvrtsys.ms.fixed']
for vis in datawith1924:
for vis in datawith1924:
   setjy(vis=vis,field='1924-292',fluxdensity=flux1924,usescratch=F)
   setjy(vis=vis,field='1924-292',fluxdensity=flux1924,usescratch=F)
Line 835: Line 1,008:
   setjy(vis=vis,field='nrao530*',fluxdensity=fluxnrao530,usescratch=F)
   setjy(vis=vis,field='nrao530*',fluxdensity=fluxnrao530,usescratch=F)


setjy(vis='uid___A002_X3d55cb_X575.antwvrtsys.ms',field='3c279',
datawith3c279=['uid___A002_X3d55cb_X575.antwvrtsys.ms.fixed']
      fluxdensity=flux3c279,usescratch=F)
for vis in datawith3c279:
setjy(vis='uid___A002_X3d55cb_X575.antwvrtsys.ms',field='1625-254',
  setjy(vis=vis,field='3c279',fluxdensity=flux3c279,usescratch=F)
      fluxdensity=flux1625,usescratch=F)
  setjy(vis=vis,field='1625-254',fluxdensity=flux1625,usescratch=F)
setjy(vis='uid___A002_X3d55cb_X575.antwvrtsys.ms',field='nrao530*',
  setjy(vis=vis,field='nrao530*',fluxdensity=fluxnrao530,usescratch=F)
      fluxdensity=fluxnrao530,usescratch=F)
</source>
</source>


Line 848: Line 1,020:
#This new amplitude calibration will be used in the applycal.
#This new amplitude calibration will be used in the applycal.


data=['uid___A002_X3d4118_X39b.antwvrtsys.ms',
data=['uid___A002_X3d4118_X39b.antwvrtsys.ms.fixed',
       'uid___A002_X3d55cb_X575.antwvrtsys.ms',
       'uid___A002_X3d55cb_X575.antwvrtsys.ms.fixed',
       'uid___A002_X3d55cb_X90c.antwvrtsys.ms',
       'uid___A002_X3d55cb_X90c.antwvrtsys.ms.fixed',
       'uid___A002_X3d55cb_Xb50.antwvrtsys.ms']
       'uid___A002_X3d55cb_Xb50.antwvrtsys.ms.fixed']
# Match up intents with source names
# Match up intents with source names
pcal='1625-254'
pcal='1625-254'
Line 862: Line 1,034:
           '1924-292,Juno,1625-254,nrao530*',
           '1924-292,Juno,1625-254,nrao530*',
           '1924-292,Juno,1625-254,nrao530*']
           '1924-292,Juno,1625-254,nrao530*']
calchan='0~3:20~3820'
refant='DV14'


for vis in range(len(data)):
for vis in range(len(data)):
Line 869: Line 1,044:
         gaintable=[data[vis]+'.bandpass_bpoly.bcal',data[vis]+'.intphase.gcal'])
         gaintable=[data[vis]+'.bandpass_bpoly.bcal',data[vis]+'.intphase.gcal'])
</source>
</source>


== Application of calibration tables ==
== Application of calibration tables ==
Line 876: Line 1,050:


<source lang="python">
<source lang="python">
data=['uid___A002_X3d4118_X39b.antwvrtsys.ms',
data=['uid___A002_X3d4118_X39b.antwvrtsys.ms.fixed',
       'uid___A002_X3d55cb_X575.antwvrtsys.ms',
       'uid___A002_X3d55cb_X575.antwvrtsys.ms.fixed',
       'uid___A002_X3d55cb_X90c.antwvrtsys.ms',
       'uid___A002_X3d55cb_X90c.antwvrtsys.ms.fixed',
       'uid___A002_X3d55cb_Xb50.antwvrtsys.ms']
       'uid___A002_X3d55cb_Xb50.antwvrtsys.ms.fixed']
# Match up intents with source names
# Match up intents with source names
pcal='1625-254'
pcal='1625-254'
Line 928: Line 1,102:


== Plot corrected data ==
== Plot corrected data ==
[[File:uid___A002_X3d55cb_X90c.antwvrtsys.ms.fixed.cal.time.spw1.png|200px|thumb|right|'''Fig. 11.''' This plot shows amplitude (flux) versus time for the dataset X90c in spw 1.]]
[[File:uid___A002_X3d55cb_X90c.antwvrtsys.ms.fixed.cal.time.phase.spw1.png|200px|thumb|right|'''Fig. 12.''' Phase vs time for all the sources in X90c (spw 1).]]
[[File:uid___A002_X3d4118_X39b.antwvrtsys.ms.fixed.cal.freq.amp.1625-254.png|200px|thumb|right|'''Fig. 13.''' Amplitude (flux) vs frequency for 1625-254 in X39b.]]
[[File:uid___A002_X3d55cb_X575.antwvrtsys.ms.fixed.cal.freq.phase.3c279.png|200px|thumb|right|'''Fig. 14.''' Phase vs frequency for 3c279 in X575.]]
[[File:uid___A002_X3d4118_X39b.antwvrtsys.ms.fixed.science.freq.amp1.png|200px|thumb|right|'''Fig. 15.''' Amplitude vs frequency for the science target, colored by fields for X39b and spw 1, showing CO (6-5).]]
[[File:uid___A002_X3d4118_X39b.antwvrtsys.ms.fixed.science.freq.amp3.png|200px|thumb|right|'''Fig. 16.''' Amplitude vs frequency as for Figure 15, but for spw 3.]]


The next commands will help you visualize the result of the application of the calibration tables to the data. You can check if the amplitudes and phases vs time and frequency look reasonable for all the sources, in particular for the science target fields.
The next commands will help you visualize the result of the application of the calibration tables to the data. You can check if the amplitudes and phases vs time and frequency look reasonable for all the sources, in particular for the science target fields.


The next command will produce four plots, as the one we showed in Figure 4, but the amplitude in this new plots corresponds to flux because we now have calibrated data.
The next command will produce four plots, as the one we showed in Figure 4, but the amplitude in this new plots corresponds to flux because we now have calibrated data.
See Figure 10 for an example of it. It is important to check that all the sources have similar amplitude (flux) in the different spws and datasets.
See Figure 11 for an example of it. It is important to check that all the sources have similar amplitude (flux) in the different spws and datasets.
 
[[File:uid___A002_X3d55cb_X90c.antwvrtsys.ms.cal.time.spw1.png|200px|thumb|right|'''Fig. 10.''' This plot shows amplitude (flux) versus time for the dataset X90c in spw 1.]]


<source lang="python">
<source lang="python">
data=['uid___A002_X3d4118_X39b.antwvrtsys.ms',
data=['uid___A002_X3d4118_X39b.antwvrtsys.ms.fixed',
       'uid___A002_X3d55cb_X575.antwvrtsys.ms',
       'uid___A002_X3d55cb_X575.antwvrtsys.ms.fixed',
       'uid___A002_X3d55cb_X90c.antwvrtsys.ms',
       'uid___A002_X3d55cb_X90c.antwvrtsys.ms.fixed',
       'uid___A002_X3d55cb_Xb50.antwvrtsys.ms']
       'uid___A002_X3d55cb_Xb50.antwvrtsys.ms.fixed']
# Match up intents with source names
# Match up intents with source names
pcal='1625-254'
pcal='1625-254'
Line 966: Line 1,150:


It is also important to check the phases vs time for all the sources. The next command will get you the corresponding plots, four for each dataset. You can see that while the bandpass and amplitude calibrator have very concentrated phases around 0 degrees, the phase calibrator and the science target do not.
It is also important to check the phases vs time for all the sources. The next command will get you the corresponding plots, four for each dataset. You can see that while the bandpass and amplitude calibrator have very concentrated phases around 0 degrees, the phase calibrator and the science target do not.
In Figure 11 we show an example of these plots for spw 3 in the dataset X90c.
In Figure 12 we show an example of these plots for spw 3 in the dataset X90c.
 
[[File:uid___A002_X3d55cb_X90c.antwvrtsys.ms.cal.time.phase.spw1.png|200px|thumb|right|'''Fig. 11.''' Phase vs time for all the sources in X90c (spw 1).]]


<source lang="python">
<source lang="python">
Line 979: Line 1,161:
</source>
</source>


We now check the amplitude of the sources vs frequency. This is important since we expect that all the spw have very similar behavior. You will have one plot for each field for each dataset. In Figure 12 we show the case for 1625-254, our phase calibrator, for X90c.
We now check the amplitude of the sources vs frequency. This is important since we expect that all the spw have very similar behavior. You will have one plot for each field for each dataset. In Figure 13 we show the case for 1625-254, our phase calibrator, for X90c.
 
[[File:uid___A002_X3d4118_X39b.antwvrtsys.ms.cal.freq.amp.1625-254.png|200px|thumb|right|'''Fig. 12.''' Amplitude (flux) vs frequency for 1625-254 in X39b.]]


<source lang="python">
<source lang="python">
Line 992: Line 1,172:
       plotfile='aftercal_plots/'+data[vis]+'.cal.freq.amp.'+field+'.png')
       plotfile='aftercal_plots/'+data[vis]+'.cal.freq.amp.'+field+'.png')
</source>
</source>
The next command will produce similar plots but this time of phase vs frequency (see Figure 13 for an example).
The next command will produce similar plots but this time of phase vs frequency (see Figure 14 for an example).
You will notice that only strong sources, like 3c279, will show clearly phases concentrated around 0 degrees.
You will notice that only strong sources, like 3c279, will show clearly phases concentrated around 0 degrees.
[[File:uid___A002_X3d55cb_X575.antwvrtsys.ms.cal.freq.phase.3c279.png|200px|thumb|right|'''Fig. 13.''' Phase vs frequency for 3c279 in X575.]]


<source lang="python">
<source lang="python">
Line 1,006: Line 1,184:
       plotfile='aftercal_plots/'+data[vis]+'.cal.freq.phase.'+field+'.png')
       plotfile='aftercal_plots/'+data[vis]+'.cal.freq.phase.'+field+'.png')
</source>
</source>
Finally, for our science target, we plot amplitude (flux) vs frequency for all the spectral windows. See an example of this is Figure 14 and 15.


[[File:IRAS16293B9_spw3_X90c_cal.png|200px|thumb|right|'''Fig. 14.''' Amplitude vs frequency for the science target, colored by fields. This corresponds to X90c in the spw 3.]]
Finally, for our science target, we plot amplitude (flux) vs frequency for all the spectral windows. See an example of this is Figure 15 and 16. If you look carefully at these plots, you will notice that the datasets at low elevation show much more line emission - this is due to the  
 
shorter projected baselines present for these datasets (X39b and Xb50). In contrast, you start to see some weak absorption for the dataset at high elevation: X90c.  
[[File:IRAS16293B9_spw1_X90c_cal.png|200px|thumb|right|'''Fig. 15.''' Amplitude vs frequency as for Figure 14, but for spw 1 showing the CO (6-5) emission line.]]


<source lang="python">
<source lang="python">
os.system('rm -rf aftercal_plots/*science.freq.amp*.png')
os.system('rm -rf aftercal_plots/*science.freq.amp*.png')
for vis in range(len(data)):
for vis in range(len(data)):
  for spw in SPW:   
     plotms(vis=data[vis],field=science,xaxis='freq', yaxis='amp',
     plotms(vis=data[vis],field=science,xaxis='freq', yaxis='amp',
       spw=spw,avgtime='1e8',avgscan=T,
       spw=spw,avgtime='1e8',avgscan=T,
       coloraxis='field',xselfscale=T,ydatacolumn='corrected',
       coloraxis='field',xselfscale=T,ydatacolumn='corrected',
       plotfile='aftercal_plots/'+data[vis]+'.science.freq.amp.png')  
       plotfile='aftercal_plots/'+data[vis]+'.science.freq.amp'+spw+'.png',
      title=data[vis].split('_')[-1].split('.')[0]+'.IRAS16293_spw'+spw)
</source>
</source>


Line 1,039: Line 1,217:
       coloraxis='field',xselfscale=T,ydatacolumn='corrected')
       coloraxis='field',xselfscale=T,ydatacolumn='corrected')
</source>
</source>


== Split and concatenate the calibrated data ==
== Split and concatenate the calibrated data ==
Line 1,047: Line 1,224:
<source lang="python">
<source lang="python">
# Splitting final calibrated datasets
# Splitting final calibrated datasets
data=['uid___A002_X3d4118_X39b.antwvrtsys.ms',
data=['uid___A002_X3d4118_X39b.antwvrtsys.ms.fixed',
       'uid___A002_X3d55cb_X575.antwvrtsys.ms',
       'uid___A002_X3d55cb_X575.antwvrtsys.ms.fixed',
       'uid___A002_X3d55cb_X90c.antwvrtsys.ms',
       'uid___A002_X3d55cb_X90c.antwvrtsys.ms.fixed',
       'uid___A002_X3d55cb_Xb50.antwvrtsys.ms']
       'uid___A002_X3d55cb_Xb50.antwvrtsys.ms.fixed']


for vis in data:
for vis in data:
   split(vis=vis,outputvis='%s.cal.IRAS16293.ms'%(vis.split('.')[0]),
   split(vis=vis,outputvis='%s.cal.IRAS16293.ms.fixed'%(vis.split('.')[0]),
     datacolumn='corrected',field='IRAS16293*',keepflags=False)
     datacolumn='corrected',field='IRAS16293*',keepflags=False)
</source>
</source>
Line 1,061: Line 1,238:
<source lang="python">
<source lang="python">
# Concatenating the final split files
# Concatenating the final split files
data=['uid___A002_X3d4118_X39b.cal.IRAS16293.ms',
concatdata=['uid___A002_X3d4118_X39b.cal.IRAS16293.ms.fixed',
     'uid___A002_X3d55cb_X575.cal.IRAS16293.ms',
     'uid___A002_X3d55cb_X575.cal.IRAS16293.ms.fixed',
     'uid___A002_X3d55cb_X90c.cal.IRAS16293.ms',
     'uid___A002_X3d55cb_X90c.cal.IRAS16293.ms.fixed',
     'uid___A002_X3d55cb_Xb50.cal.IRAS16293.ms']
     'uid___A002_X3d55cb_Xb50.cal.IRAS16293.ms.fixed']


concat(vis=data,concatvis='IRAS16293_Band9.ms')
concat(vis=concatdata,concatvis='IRAS16293_Band9.fixed.ms')
</source>
 
To speed up imaging we apply 60 second time averaging to the concatenated dataset.
 
<source lang="python">
# 60s time averaging
split(vis='IRAS16293_Band9.fixed.ms', datacolumn='data', timebin='60s',
      outputvis='IRAS16293_Band9.fixed.rebin.ms')
</source>
 
As a final step, we zero the rows of the pointing table because it is
quite large and is not currently needed by the imaging software for
mosaics -indeed its presence will cause an error during imaging if you
skip this step.
 
<source lang="python">
# Remove rows of pointing table
tb.open('IRAS16293_Band9.fixed.rebin.ms/POINTING', nomodify = False)
a = tb.rownumbers()
tb.removerows(a)
tb.close()  
</source>
</source>


Now you have everything you need to the imaging part. Follow [http://casaguides.nrao.edu/index.php?title=IRAS16293_Band9_-_Imaging '''IRAS16293 Band9 - Imaging'''] to go to the imaging section of this casaguide.
Now you have completed the calibration and have everything you need to carry out the imaging stage. Follow [http://casaguides.nrao.edu/index.php?title=IRAS16293_Band9_-_Imaging_for_CASA_3.4 '''IRAS16293 Band9 - Imaging'''] to go to the imaging section of this casaguide.


{{Checked 3.4.0}}
{{Checked 3.4.0}}

Latest revision as of 17:49, 15 June 2012


WARNING: On June 15, 2012 the calibration guide and the final data products (calibrated science 
data: IRAS16293_Band9_CalibratedMS_FIXED.tgz and reference images: IRAS16293_Band9_ReferenceImages_FIXED.tgz)) 
were changed to correct for a 1.2" position error in the phase calibrator (1625-254). Without 
correction, the science images will suffer from a similar offset. 

Overview

This part of the casa guide will guide you through the basic inspection of the data, paying special attention to identifying data that needs to be flagged. The guide shows the complete process to get fully calibrated data.

The general procedure in this guide follows the other ALMA CASA guides: NGC3256Band3 and AntennaeBand7.


Unpack the Data

Once you have downloaded the IRAS16293_Band9_UnCalibratedMS.tgz, unpack the file in a terminal outside CASA using

tar -xvzf IRAS16293_Band9_UnCalibratedMS.tgz

cd IRAS16293_Band9_UnCalibratedMS

You have a number of files with extensions ".ms", which are CASA measurement set (MS) files. You will also see files containing system temperature (Tsys), water vapor radiometer (WVR), and antenna position information.

To start CASA type

casapy

Be sure that you are using the right version indicated for this guide.

Confirm your version of CASA

This guide has been written for CASA release 3.4.0. Please confirm your version before proceeding. This guide will not work in earlier versions (< r19874).

# In CASA
version = casalog.version()
print "You are using " + version
if (int(version.split()[4][1:-1]) < 19874):
    print "\033[91m YOUR VERSION OF CASA IS TOO OLD FOR THIS GUIDE."
    print "\033[91m PLEASE UPDATE IT BEFORE PROCEEDING."
else:
    print "Your version of CASA is appropriate for this guide."

Install Analysis Utilities

Analysis Utilities (or analysisUtils for short) is a small set of Python scripts that provide a number of analysis and plotting utilities for ALMA data reduction. This guide uses a few of these utilities. They are very easy to install (just download and untar). See

http://casaguides.nrao.edu/index.php?title=Analysis_Utilities

for a full description and download instructions. If you do not wish to do this, see a CASA 3.3 version of one of the other ALMA guides for alternative (but slow) plotting options. Analysis Utilities are updated frequently so if its been a while since you installed it, it's probably worth doing it again. If you are at an ALMA site or ARC, the analysis utilities are probably already installed and up to date.

Initial Inspection

The first step we will do through all the calibration process is to define an array with the uid's that correspond to the dataset names. This will allow us to make the calibration of the four datasets one after another, using a for-loop inside python. We will then calibrate the data individually and concatenate them at the end, before proceeding with the imaging part.

Note that if you exit CASA and want to continue with the calibration using these arrays, you will have to re-issue the command again to make it available for the current CASA execution.

To start, and give an example of this process, we will create txt format files for the output of the listobs task, which will give us useful information about the observations.

# Array containing the uid names

rawdata=['uid___A002_X3d4118_X39b.ms','uid___A002_X3d55cb_X575.ms',
         'uid___A002_X3d55cb_Xb50.ms','uid___A002_X3d55cb_X90c.ms']

# We create the text files for listobs for each dataset

for data in rawdata:
  listobs(vis=data,listfile=data+'.listobs')

Note that after cutting and pasting a for-loop you often have to press return several times to make the command execute. The output will be sent to the CASA logger. Next, we show an example of a useful part of the output that the first listobs command produces.

Fields: 11
  ID   Code Name                RA              Decl          Epoch   SrcId nVis   
  0    none 1924-292            19:24:51.05600 -29.14.30.1280 J2000   0     169125 
  1    none nrao530 ph          17:33:02.72400 -13.04.49.4860 J2000   1     289170 
  2    none Juno                16:25:31.63031 -05.49.08.9209 J2000   2     82890  
  3    none 1625-254            16:25:46.98000 -25.27.38.3300 J2000   3     276480 
  4    none IRAS16293-2422-a    16:32:22.99200 -24.28.36.0000 J2000   4     132450 
  5    none IRAS16293-2422-a    16:32:22.47925 -24.28.36.0000 J2000   4     99915  
  6    none IRAS16293-2422-a    16:32:22.73563 -24.28.36.0000 J2000   4     99960  
  7    none IRAS16293-2422-a    16:32:22.73563 -24.28.32.5000 J2000   4     99915  
  8    none IRAS16293-2422-a    16:32:22.47925 -24.28.29.0000 J2000   4     99945  
  9    none IRAS16293-2422-a    16:32:22.73563 -24.28.29.0000 J2000   4     99945  
  10   none IRAS16293-2422-a    16:32:22.99200 -24.28.29.0000 J2000   4     99915  
   (nVis = Total number of time/baseline visibilities per field) 
Spectral Windows:  (25 unique spectral windows and 2 unique polarization setups)
  SpwID  #Chans Frame Ch1(MHz)    ChanWid(kHz)  TotBW(kHz)  Corrs  
  0           4 TOPO  184550      1500000       7500000     I   
  1         128 TOPO  231257.813  15625         2000000     XX  YY  
  2           1 TOPO  232234.375  1796875       1796875     XX  YY  
  3         128 TOPO  229257.813  15625         2000000     XX  YY  
  4           1 TOPO  230234.375  1796875       1796875     XX  YY  
  5         128 TOPO  217242.188  15625         2000000     XX  YY  
  6           1 TOPO  216234.375  1796875       1796875     XX  YY  
  7         128 TOPO  215242.188  15625         2000000     XX  YY  
  8           1 TOPO  214234.375  1796875       1796875     XX  YY  
  9         128 TOPO  703257.813  15625         2000000     XX  YY  
  10          1 TOPO  704234.375  1796875       1796875     XX  YY  
  11        128 TOPO  692492.188  15625         2000000     XX  YY  
  12          1 TOPO  691484.375  1796875       1796875     XX  YY  
  13        128 TOPO  690492.188  15625         2000000     XX  YY  
  14          1 TOPO  689484.375  1796875       1796875     XX  YY  
  15        128 TOPO  688492.188  15625         2000000     XX  YY  
  16          1 TOPO  687484.375  1796875       1796875     XX  YY  
  17       3840 TOPO  703312.744  488.28125     1875000     XX  YY  
  18          1 TOPO  704249.756  1875000       1875000     XX  YY  
  19       3840 TOPO  692237.256  488.28125     1875000     XX  YY  
  20          1 TOPO  691299.756  1875000       1875000     XX  YY  
  21       3840 TOPO  690437.256  488.28125     1875000     XX  YY  
  22          1 TOPO  689499.756  1875000       1875000     XX  YY  
  23       3840 TOPO  688437.256  488.28125     1875000     XX  YY  
  24          1 TOPO  687499.756  1875000       1875000     XX  YY  
  
Antennas: 15:
  ID   Name  Station   Diam.    Long.         Lat.         
  0    DA41  A003      12.0 m   -067.45.16.5  -22.53.27.0  
  1    DA43  A075      12.0 m   -067.45.17.9  -22.53.21.4  
  2    DA47  A026      12.0 m   -067.45.18.8  -22.53.28.3  
  3    DV02  A077      12.0 m   -067.45.10.1  -22.53.25.9  
  4    DV03  A137      12.0 m   -067.45.15.2  -22.53.22.7  
  5    DV05  A082      12.0 m   -067.45.08.3  -22.53.29.2  
  6    DV07  A076      12.0 m   -067.45.20.5  -22.53.33.8  
  7    DV09  A046      12.0 m   -067.45.17.0  -22.53.29.3  
  8    DV10  A071      12.0 m   -067.45.19.9  -22.53.23.5  
  9    DV12  A011      12.0 m   -067.45.14.4  -22.53.28.4  
  10   DV13  A072      12.0 m   -067.45.12.6  -22.53.24.0  
  11   DV14  A025      12.0 m   -067.45.18.7  -22.53.27.4  
  12   DV15  A074      12.0 m   -067.45.12.1  -22.53.32.0  
  13   DV17  A138      12.0 m   -067.45.17.1  -22.53.34.4  
  14   DV18  A053      12.0 m   -067.45.17.3  -22.53.31.2  

In the previous output you can see the ID that is assigned to each source, starting with the number 0. 1924-292 and 3c279 are the calibrators for bandpass, Juno for amplitude (flux), 1625-254 is our phase calibrator and nrao530 ph serves as a check (phase calibrator) source. The remaining 7 fields of IRAS16293-2422-a are the seven pointings for our mosaic of the target source.

Spectral windows are also marked with numbers from 0 to 24, with number 0 containing WVR information. Spws 17, 19, 21, and 23 contain the science data (FDM mode). The CO (6-5) line emission is contained in spw 19. Spw 18, 20, 22, and 24 contain channel averages of the data from spectral windows 17, 19, 21, 23, respectively. These spws will not be used for the offline data reduction. All the remaining spws that appear in the section of "Sources", and that do not appear in the "Spectral Windows" section are related to WVR measurements for each antenna, so you will not need them for the calibration either. Spws 9, 11, 13, and 15 are associated with tsys measurements, and we will apply these measurements to the science spws later.

Finally before we go further we explicitly save the current flag state of the data. If you ever decide to start completely over, you should substitute 'restore' for 'save' in the command below to restore the flag state to its original value.

# Use flagmanager to save current flag state.

for vis in rawdata:
  flagmanager(vis=vis,mode='save',versionname='OriginalFlagState')

Generation and visualization of the antenna positions, Tsys and WVR tables

Next we will generate the tables we need to apply to the data: antenna position, system temperature and water vapor radiometer. Once the tables are generated we will produce plots of them and inspect them to make sure whether they have issues that might affect their application to the data. Whenever we see an odd behavior in the tables we need to flag the corresponding science data to prevent wrong results in the calibration steps.

System Temperature

We produce the Tsys tables with the next command

# Create Tsys
os.system('rm -rf *tdm.tsys')
for vis in rawdata:
    print "Creating TDM Tsys Table for "+vis
    gencal(vis=vis,caltable=vis+'.tdm.tsys',spw='9,11,13,15',caltype='tsys')

The next command, that comes from the Analysis Utils package will plot the tsys in the next way: it will produce many plots, each one of them will show an antenna, with the four spws that the tsys tables cover, for all the targets, and with different colors for different times, so you can trace the behavior for tsys with time, among others. Note that in spw 19, the overlap with the tsys spw (11) is not set correctly. This is due to an error in the frequencies for the tsys when the observations were done. You do not have to worry about this, since any issue coming from that error has already been fixed. Note, however that the portions of the spectra that do not have tsys information cannot be used. This does not represent a problem, since that part corresponds to the edge of the baseband. Also note that the CO (6-5) line is not affected by this. In Figure 1 you will see the corresponding plot for one of the datasets (X90c) showing antenna 0 (DA41).

Fig. 1. Example for the output of the command that plots the tsys spws. Find a description in the text.


# Plot TDM Tsys tables with times overlayed for each antenna
os.system('rm -rf Tsysplots/*time*')
for vis in rawdata:
  aU.plotbandpass(caltable=vis+'.tdm.tsys',ms=vis,
                overlay='time', xaxis='freq', showatm=True,
                yaxis='amp',subplot=22,interactive=False,
                chanrange='5~122',showfdm=True,
                figfile='Tsysplots/'+vis+'.time.tdm.tsys.png')

Another useful set of plots is generated by the next command. Each resulting plot shows all the antennas in a scan, for one spectral windows. This is useful to check Tsys for all the antennas in a certain time. In Figure 2 you can see an example.

Fig. 2. Second sample plot for tsys plotting. The first four scans of the observations for spw 17 are shown, with all the antennas per plot.
# Plot TDM Tsys tables with antennas overlayed for each time 
os.system('rm -rf Tsysplots/*antenna*')
for vis in rawdata:
  aU.plotbandpass(caltable=vis+'.tdm.tsys',ms=vis,
                overlay='antenna', xaxis='freq', showatm=True,
                yaxis='amp',subplot=22,interactive=False,
                chanrange='5~122',showfdm=True,
                figfile='Tsysplots/'+vis+'.antenna.tdm.tsys.png')

Go through all the plots and make sure you notice all the next issues, since we will need to flag the corresponding science data.

 X90c
 DV05 ripples all spw  
 Otherwise 600 to 1200 47 to 57% transmission

 X575
 Otherwise 1300 to 3000 28 to 39% transmission
 DV05 ripples all spw 

 Xb50
 800 to 2500 37 to 47% transmission
 DV05 ripples all spw 

 X39b
 500 to 800 56 to 65%
 DA43 Tsys extremely high for spw=23
 DV05 ripples all spw and one time with bad YY 
 DV18 extremely high for spw=23

Water Vapor Radiometer

We now generate the WVR tables by executing

rawdata=['uid___A002_X3d4118_X39b.ms','uid___A002_X3d55cb_X575.ms',
         'uid___A002_X3d55cb_Xb50.ms','uid___A002_X3d55cb_X90c.ms']

# Other values can be left at defaults
os.system('rm -rf *.wvrgcal')
for vis in rawdata:
    wvrgcal(vis=vis,caltable=vis+'.wvr',toffset=-1)

The integration time for the science spectral windows is 6.048 seconds. So we need to average the WVR 1 second solutions to match.

# Smoothing the data
for vis in rawdata:
  print "Smoothing wvr table for  "+vis
  smoothcal(vis=vis,tablein=vis+'.wvr',caltable=vis+'.wvr.smooth',
            smoothtype='mean',smoothtime=6.048)

Now, for the plotting of the wvr tables, we will employ again the analysis utils. This command will create a directory with all the plots inside, each one of them corresponding to different datasets, baselines and targets, using different colors. In Figure 3 you can see an example of the output for spectral window 1. Note that the command below only creates the plots for that spw, since the others are the same except for a scale factor that is the ratio of frequencies. Note that in all datasets, DV15 has bad wvr behavior.

Fig. 3. Phase corrections as function of time for the dataset X90c, where you can see an example of the odd behavior of DV15-related baselines.
# Plotting wvr tables for spw 1
os.system('rm -rf WVRplots')
for vis in rawdata:
  aU.plotWVRSolutions(caltable=vis+'.wvr.smooth',
              yrange=[-180,180],figfile='WVRplots/'+vis+'.wvr.smooth.png',
              ms=vis,spw='1',interactive=False)

Antenna position

Finally, we need to create the tables to correct for the antenna positions. For that, execute

# Generate the tables for antenna positions
rawdata=['uid___A002_X3d4118_X39b.ms','uid___A002_X3d55cb_X575.ms',
         'uid___A002_X3d55cb_Xb50.ms','uid___A002_X3d55cb_X90c.ms']

# For the data taken on April 16
antennas='DV18,DV12,DV09,DV10,DV13,DV05,DA41,DV14,DA43,DV17,DV07,DV15'
parameter=[0.00684754783288,-0.0188271608204,-0.00732900947332,
0.000481227821554,-2.48519197708e-05,-0.00048066949752,
0.000741523189806,-0.000843518709779,-0.000409113120338,
0.00111751207514,-0.000891124202747,-0.000492116263895,
0.00100313175518,-0.00084909900337,-0.000594194887412,
0.00105206234994,-0.000857719516042,-0.00059565765641,
0.00060776527971,-0.000525096431375,-0.000386487226933,
0.000521688542419,-0.000144920371688,-0.000258341151538,
0.00104740676978,-0.000899314162014,-0.000483014592547,
0.000920921096116,-0.00106288533016,-0.000546355296173,
0.000867350985729,-0.000990079062422,-0.00038123971161,
0.000567245762795,-0.000422531738877,-0.000741700641811]
gencal(vis = 'uid___A002_X3d4118_X39b.ms',
  caltable = 'uid___A002_X3d4118_X39b.ms.antpos',
  caltype = 'antpos',antenna=antennas,parameter=parameter)

# and for April 17
antennas='DV05,DV09,DV10,DV13,DV12,DA41,DV14,DA43,DV17,DV07,DV15'
parameter=[0.00105206234994,-0.000857719516042,-0.00059565765641,
0.000741523189806,-0.000843518709779,-0.000409113120338,
0.00111751207514,-0.000891124202747,-0.000492116263895,
0.00100313175518,-0.00084909900337,-0.000594194887412,
0.000481227821554,-2.48519197708e-05,-0.00048066949752,
0.00060776527971,-0.000525096431375,-0.000386487226933,
0.000521688542419,-0.000144920371688,-0.000258341151538,
0.00104740676978,-0.000899314162014,-0.000483014592547,
0.000920921096116,-0.00106288533016,-0.000546355296173,
0.000867350985729,-0.000990079062422,-0.00038123971161,
0.000567245762795,-0.000422531738877,-0.000741700641811]
for vis in range(1,4):
  gencal(vis=rawdata[vis],
    caltable=rawdata[vis]+'.antpos',
    caltype ='antpos',antenna=antennas,parameter=parameter)

Now, based on the behavior of the tsys and wvr tables, we will flag the corresponding data, using the next commands. You can employ similar executions to flag other data you might want to remove.

rawdata=['uid___A002_X3d4118_X39b.ms','uid___A002_X3d55cb_X575.ms',
         'uid___A002_X3d55cb_Xb50.ms','uid___A002_X3d55cb_X90c.ms']


# Flagging corresponding science data for tsys and wvr showing problems
for vis in rawdata:
  flagdata(vis=vis,autocorr = T,flagbackup = F)
  flagdata(vis=vis,mode='shadow',flagbackup=F)
  flagdata(vis=vis,antenna='DV05,DV15',flagbackup=F)

vis='uid___A002_X3d4118_X39b.ms'
flagdata(vis=vis,antenna='DA43,DV18',spw='23',flagbackup=F)

Applying antpos, tsys, and wvr tables and splitting the data

To apply the Tsys tables we need to separate sources that have Tsys measurements of their own from those that do not. Since only one IRAS16293 field (id=4) has Tsys, it goes into the "not" category. As you probably noted from the messages from plotbandpass, the Tsys measurements on Juno did not yield usable values, causing them to be flagged, so we will apply the nearest source in elevation to it: IRAS16293 field id=4. This can be checked with the plots generated by the following plotms commands. In Figure 4 you can see the output of this command for Xb50.

Fig. 4. Elevation for all the sources versus time for Xb50. Note that this dataset has very low elevations.
# Re-entering our array
rawdata=['uid___A002_X3d4118_X39b.ms','uid___A002_X3d55cb_X575.ms',
         'uid___A002_X3d55cb_Xb50.ms','uid___A002_X3d55cb_X90c.ms']

# Elevation plots to understand what the elevation range for each dataset is.
for vis in rawdata:
  plotms(vis=vis, 
       field='',xaxis='time', yaxis='elevation',antenna='',
       spw='17', avgchannel='3840',coloraxis='field',
       ydatacolumn='data',plotfile=vis+'elevation.png',title=vis)

To apply the Tsys we need to explicitly tell applycal which TDM (128 channels) spws go with which FDM spw (3840 channels). An entry is needed for every spw, because the position in the list corresponds to the spw id to be corrected. There is a helper function to assist you in figuring this out.

from recipes.almahelpers import tsysspwmap
tsysspwmap(vis='uid___A002_X3d4118_X39b.ms',
           tsystable='uid___A002_X3d4118_X39b.ms.tdm.tsys')

# This will print:
# [0,1,2,3,4,5,6,7,8,9,9,11,11,13,13,15,15,9,9,19,20,13,13,15,15]

The helper gets a little confused in this case with spws 19 and 20 because the FDM was not perfectly positioned within the TDM window as you can see from the Tsys plots. In a case like this, we need to put in by hand the TDM window that should match the FDM window.

As you could see from the first run of listobs, 1924-292 is the bandpass calibrator for three of the datasets, and 3c279 is for one of them. For this reason, the application of the tables is split into two parts.

For the applycal task, note that the first value in the interp parameter for each gaintable gives the desired time interpolation type. The second parameter indicates the desired frequency axis interpolation. If the second value is not set it is assumed to be linear, but only if the input table has a frequency axis (like Tsys and Bandpass). Spline seems to work best for the Tsys TDM to FDM frequency interpolation.

rawdata=['uid___A002_X3d4118_X39b.ms','uid___A002_X3d55cb_X575.ms',
         'uid___A002_X3d55cb_Xb50.ms','uid___A002_X3d55cb_X90c.ms']

tsysspw=[0,1,2,3,4,5,6,7,8,9,9,11,11,13,13,15,15,9,9,11,11,13,13,15,15]

# Datasets with 1924-292 as the bandpass calibrator
for vis in [rawdata[0],rawdata[2],rawdata[3]]:
  print "Applying calibration for "+vis
  field_Tsys=['1924-292','nrao530 ph']
  for field in field_Tsys:
    print "For Field= "+field
    applycal(vis=vis,field=field,
             spw = '17,19,21,23',
             gaintable = [vis+'.tdm.tsys', 
                          vis+'.wvr.smooth', 
                          vis+'.antpos'],
             gainfield = [field,field,''],
             spwmap=[tsysspw,[],[]],
             interp = ['linear,spline','nearest',''],calwt = T,
             flagbackup = F)


# This next applycal takes care of the fact that one dataset has 3C279
# instead of 1924-292 as the bandpass calibrator 

for vis in [rawdata[1]]:
  print "Applying calibration for "+vis
  field_Tsys=['3c279','nrao530 ph']
  for field in field_Tsys:
    print "For Field= "+field
    applycal(vis=vis,field=field,
             spw = '17,19,21,23',
             gaintable = [vis+'.tdm.tsys', 
                          vis+'.wvr.smooth', 
                          vis+'.antpos'],
             gainfield = [field,field,''],
             spwmap=[tsysspw,[],[]],
             interp = ['linear,spline','nearest',''],calwt = T,
             flagbackup = F)

# Next we do the "noTsys" sources for all datasets

for vis in rawdata:
  print "Applying calibration for "+vis
  field_noTsys=['Juno','1625-254','IRAS16293*']
  for field in field_noTsys:
    print "For Field= "+field
    applycal(vis =vis,field=field,
             spw = '17,19,21,23',
             gaintable = [vis+'.tdm.tsys', 
                          vis+'.wvr.smooth', 
                          vis+'.antpos'],
             gainfield = ['4',field,''],
             spwmap=[tsysspw,[],[]],
             interp = ['linear,spline','nearest',''],calwt = T,
             flagbackup = F)

# Splitting the science spws

for vis in rawdata:
  split(vis=vis,outputvis=('%s.antwvrtsys.ms'%(vis.split('.')[0])),
        datacolumn='corrected',spw='17,19,21,23',keepflags=False)

Fix Phase Calibrator Position

The position used for the phase calibrator (1625-254) in these observations is offset by about 1.2" toward postive R.A from its correct position. Below we correct the data for this offset. It is best to do this step before doing any of the gain calibration. It is notable that simply shifting the images for this offset will not give as correct a result due to phase delays caused by sky rotation during the course of each observation.

data=['uid___A002_X3d4118_X39b.antwvrtsys.ms',
      'uid___A002_X3d55cb_X575.antwvrtsys.ms',
      'uid___A002_X3d55cb_X90c.antwvrtsys.ms',
      'uid___A002_X3d55cb_Xb50.antwvrtsys.ms']

List the data to see current position of 1625-254

for vis in data:
  listobs(vis=vis,listfile=vis+'.listobs',verbose=True)

Incorrect position from listobs: 1625-254 16h25m46.98000s -25d27m38.3300s J2000

From the [http://www.vla.nrao.edu/astro/calib/manual/csource.html EVLA calibrator manual] we find that the correct position is: 16h25m46.891639s -25d27m38.326880s J2000

The task fixvis still requires scratch columns, so we must create these first (this takes a while).

for vis in data:
  clearcal(vis=vis)

Next use the task fixvis to correct the position of the calibrator in the data and the header. It will also recalculate UVWs, but this correction is very small for a 1.2" shift.

for vis in data:
  fixvis(vis=vis,outputvis=vis+'.fixed',field='1625-254',reuse=False,
         phasecenter='J2000 16h25m46.891639s -25d27m38.326880s')

Define the new dataset names that we will use here on out in the calibration.

data=['uid___A002_X3d4118_X39b.antwvrtsys.ms.fixed',
      'uid___A002_X3d55cb_X575.antwvrtsys.ms.fixed',
      'uid___A002_X3d55cb_X90c.antwvrtsys.ms.fixed',
      'uid___A002_X3d55cb_Xb50.antwvrtsys.ms.fixed']

Next, we need to remove the model data column, as fixvis currently rotates both the data AND the model data column. If you do not do the following step, you will double correct for the position offset.

for vis in data:
  delmod(vis=vis,scr=True,otf=False)

Check the new position in listobs

for vis in data:
  listobs(vis=vis,listfile=vis+'listobs_afterfix',verbose=F)

Data inspection

We now need to check for any bad behavior in the data through several plots. Once problems are identified, data can be flagged. But before that, we need to run again listobs to check that the split worked as expected. We will define our new array of split datasets, along with a list of intents that will be useful in the next steps.


# New array of datasets 
data=['uid___A002_X3d4118_X39b.antwvrtsys.ms.fixed',
      'uid___A002_X3d55cb_X575.antwvrtsys.ms.fixed',
      'uid___A002_X3d55cb_X90c.antwvrtsys.ms.fixed',
      'uid___A002_X3d55cb_Xb50.antwvrtsys.ms.fixed']

# Match up intents with source names, for the lists below, it is important 
# that the order match that of the data parameter
pcal='1625-254'
fluxcal='Juno' 
science='IRAS16293*'
check='nrao530*'
bpcal=['1924-292','3c279','1924-292','1924-292']
calfields=['1924-292,Juno,1625-254,nrao530*',
           '3c279,Juno,1625-254,nrao530*',
           '1924-292,Juno,1625-254,nrao530*',
           '1924-292,Juno,1625-254,nrao530*']

for vis in data:
    listobs(vis=vis,listfile=vis+'.listobs',verbose=True)

You can explore any of the output files by doing cat file.listobs or using any other text reader from a terminal not running CASA. Next you can see the output for X90c, and you will see the change in the spw naming, among others.

Fields: 11
  ID   Code Name                RA              Decl          Epoch   SrcId nVis   
  0    none 1924-292            19:24:51.05600 -29.14.30.1280 J2000   0     13200  
  1    none nrao530 ph          17:33:02.72400 -13.04.49.4860 J2000   1     23760  
  2    none Juno                16:25:05.61170 -05.43.27.9210 J2000   2     10560  
  3    none 1625-254            16:25:46.98000 -25.27.38.3300 J2000   3     26400  
  4    none IRAS16293-2422-a    16:32:22.99200 -24.28.36.0000 J2000   4     11880  
  5    none IRAS16293-2422-a    16:32:22.47925 -24.28.36.0000 J2000   4     11880  
  6    none IRAS16293-2422-a    16:32:22.73563 -24.28.36.0000 J2000   4     11880  
  7    none IRAS16293-2422-a    16:32:22.73563 -24.28.32.5000 J2000   4     11880  
  8    none IRAS16293-2422-a    16:32:22.47925 -24.28.29.0000 J2000   4     10560  
  9    none IRAS16293-2422-a    16:32:22.73563 -24.28.29.0000 J2000   4     10560  
  10   none IRAS16293-2422-a    16:32:22.99200 -24.28.29.0000 J2000   4     10560  
   (nVis = Total number of time/baseline visibilities per field) 
Spectral Windows:  (4 unique spectral windows and 1 unique polarization setups)
  SpwID  #Chans Frame Ch1(MHz)    ChanWid(kHz)  TotBW(kHz)  Corrs  
  0        3840 TOPO  703312.744  488.28125     1875000     XX  YY  
  1        3840 TOPO  692237.256  488.28125     1875000     XX  YY  
  2        3840 TOPO  690437.256  488.28125     1875000     XX  YY  
  3        3840 TOPO  688437.256  488.28125     1875000     XX  YY  
Sources: 20
  ID   Name                SpwId RestFreq(MHz)  SysVel(km/s) 
  0    1924-292            0     -              -            
  0    1924-292            1     -              -            
  0    1924-292            2     -              -            
  0    1924-292            3     -              -            
  1    Juno                0     -              -            
  1    Juno                1     -              -            
  1    Juno                2     -              -            
  1    Juno                3     -              -            
  2    1625-254            0     -              -            
  2    1625-254            1     -              -            
  2    1625-254            2     -              -            
  2    1625-254            3     -              -            
  3    nrao530 ph          0     -              -            
  3    nrao530 ph          1     -              -            
  3    nrao530 ph          2     -              -            
  3    nrao530 ph          3     -              -            
  4    IRAS16293-2422-a    0     -              -            
  4    IRAS16293-2422-a    1     -              -            
  4    IRAS16293-2422-a    2     -              -            
  4    IRAS16293-2422-a    3     -              -            
Antennas: 12:
  ID   Name  Station   Diam.    Long.         Lat.         
  0    DA41  A003      12.0 m   -067.45.16.5  -22.53.27.0  
  1    DA43  A075      12.0 m   -067.45.17.9  -22.53.21.4  
  2    DV02  A077      12.0 m   -067.45.10.1  -22.53.25.9  
  3    DV03  A137      12.0 m   -067.45.15.2  -22.53.22.7  
  4    DV05  A082      12.0 m   -067.45.08.3  -22.53.29.2  
  5    DV07  A076      12.0 m   -067.45.20.5  -22.53.33.8  
  6    DV09  A046      12.0 m   -067.45.17.0  -22.53.29.3  
  7    DV10  A071      12.0 m   -067.45.19.9  -22.53.23.5  
  8    DV12  A011      12.0 m   -067.45.14.4  -22.53.28.4  
  9    DV13  A072      12.0 m   -067.45.12.6  -22.53.24.0  
  10   DV14  A025      12.0 m   -067.45.18.7  -22.53.27.4  
  12   DV17  A138      12.0 m   -067.45.17.1  -22.53.34.4  

Next, we give you a set of useful plotms commands which will help you to analyze all the data in several ways. You can save a copy of the output, so you do not have to run them again every time you want to check them. This is especially useful for the ones that take a lot of time to complete. In Figure 5 you can see the output of the following plotms command, which runs with X39b. By clicking the "Next" arrow in plotms you can access the remaining spws, since the command was executed with the option iteraxis='spw'.

Fig. 5. plotms result for amplitude vs time for all the sources, which are displayed with different colors. The plot shows spw 0 for the dataset X39b.
# Check overall behavior with time
vis=data[0]
plotms(vis=vis, 
       field='',xaxis='time', yaxis='amp',antenna='',
       spw='', avgchannel='3840',coloraxis='field',
       iteraxis='spw',ydatacolumn='data',yselfscale=True)

For the next set of plotms commands, inspect each dataset, noting any problems that you notice.

# For at least one spw go antenna by antenna to look for dropouts not
# obvious in previous plot
vis=data[0]
plotms(vis=vis, 
       field='',xaxis='time', yaxis='amp',antenna='',
       spw='2', avgchannel='3840',coloraxis='field',
       iteraxis='antenna',ydatacolumn='data')

# Check out spectral properties of each source for problems
vis=data[0]
plotms(vis=vis, 
       field='',xaxis='freq', yaxis='amp',antenna='',
       spw='', avgtime='1e8',avgscan=True,coloraxis='spw',
       iteraxis='field',ydatacolumn='data',yselfscale=True)

# Examine phase of the bandpass calibrator for any problems
vis=data[0]
bp=bpcal[0] 
plotms(vis=vis, 
       field=bp,xaxis='freq', yaxis='phase',antenna='',
       spw='', avgtime='1e8',avgscan=True,avgchannel='10',coloraxis='spw',
       iteraxis='baseline',ydatacolumn='data',yselfscale=True)

Flagging

Next, based on our inspection we will proceed with the corresponding flagging. But before that, we will save the current flags state, so we can recover it later, if needed.

data=['uid___A002_X3d4118_X39b.antwvrtsys.ms.fixed',
      'uid___A002_X3d55cb_X575.antwvrtsys.ms.fixed',
      'uid___A002_X3d55cb_X90c.antwvrtsys.ms.fixed',
      'uid___A002_X3d55cb_Xb50.antwvrtsys.ms.fixed']

# Back up flag state in case you want to start over.
for vis in data:
  flagmanager(vis=vis,mode='save',versionname='Original')

# If you do start over run this first
for vis in data:
  flagmanager(vis=vis,mode='restore',versionname='Original')

Here is the flagging that we will do. If you note some other data that needs to be flagged, proceed with that as well.

# DV07 low amp after certain time. DV17 also but not as bad.
flagdata(vis='uid___A002_X3d55cb_X575.antwvrtsys.ms.fixed',
         antenna='DV07', timerange='>05:20:00',flagbackup=F)

Now, during the procedure of calibration some problems in the data will show up. Those data need to be flagged and then the calibration should be repeated. In order to save you time with this, we note here those problems, so you can flag the data now.

# PROBLEMS DISCOVERED AFTER INITIAL CALIBRATION INSPECTION

# flag low elevation scans on 1625-254 and IRAS16293
flagdata(vis='uid___A002_X3d55cb_Xb50.antwvrtsys.ms.fixed',
         timerange='>11:54:00', field='',flagbackup=F)

# flag low gains on DV02 on Juno
flagdata(vis='uid___A002_X3d55cb_Xb50.antwvrtsys.ms.fixed',
         antenna='DV02', field='Juno',flagbackup=F)

flagdata(vis='uid___A002_X3d4118_X39b.antwvrtsys.ms.fixed',
         antenna='DV02', field='Juno',flagbackup=F)

Calibration

Now we can start with the calibration itself. First, we will perform the bandpass calibration using 1924-292 and 3c279. As before, we define our list of data and match the sources with intents. Also, we will set our reference antenna (one close to the center of the array and without problems, like delays). As you can see, two different intervals for channels are used, which will make sense in a bit.

data=['uid___A002_X3d4118_X39b.antwvrtsys.ms.fixed',
      'uid___A002_X3d55cb_X575.antwvrtsys.ms.fixed',
      'uid___A002_X3d55cb_X90c.antwvrtsys.ms.fixed',
      'uid___A002_X3d55cb_Xb50.antwvrtsys.ms.fixed']
# Match up intents with source names
pcal='1625-254'
fluxcal='Juno' 
science='IRAS16293*'
check='nrao530*'
bpcal=['1924-292','3c279','1924-292','1924-292']
calfields=['1924-292,Juno,1625-254,nrao530*',
           '3c279,Juno,1625-254,nrao530*',
           '1924-292,Juno,1625-254,nrao530*',
           '1924-292,Juno,1625-254,nrao530*']

# Setup calibration parameters
prebpchan='0~3:1200~1500'
calchan='0~3:20~3820'
refant='DV14'
gaps=5
os.system('rm -rf *cal')

Note that if at some point during the calibration process you need to start over, then you will need to clear all the columns for solutions in the data, as shown next.

for vis in range(len(data)):
  delmod(vis=data[vis])
  setjy(vis=data[vis],field=fluxcal,standard='Butler-JPL-Horizons 2010',
        scalebychan=True,usescratch=False)

Now we proceed with the bandpass calibration, which will be done in two steps. In the first one, we will use a small block of channels in the center of the spws, 1200~1500, to calculate gains in phase, using a solution interval of 30 seconds. This will give us the variation of the phase throughout the observation. We will use this table later to execute the bandpass itself. Note that we are using a minimum signal to noise ratio of 2 to accept the solutions. This is low but necessary since the calibrators are not very strong in band 9.

# Calculating phase variation with time
for vis in range(len(data)):
  gaincal(vis=data[vis],caltable=data[vis]+'.bpphase.gcal',
        field=bpcal[vis],spw=prebpchan,refant=refant,
        calmode='p',solint='30s',minsnr=2.0,minblperant=4)

Next, we set our bandpass command to use our previously generated gain tables. Since we do not have a high signal to noise ratio per channel, we use a polynomial option to calculate the solutions for the bandpass. In the bandpass task, the options degamp and degphase will set the maximum degree of the polynomial that the task can use to calculate solutions. In the work log you can note what is the actual degree that the task is using. Note that the combination of solint='inf' and combine='scan' will result in one solution per scan for the calibrator.

Fig. 6. Bandpass solutions for the dataset X90c (all 4 spws) showing antenna DV09. Both amplitude and phases solutions are plotted. These data use the weaker bandpass calibrator J1924-224.
Fig. 7. Bandpass solutions for the dataset X575 (all 4 spws) showing antenna DV09. Both amplitude and phases solutions are plotted. These data use the stronger bandpass calibrator 3C279.


# Bandpass calibration using previous tables as input
for vis in range(len(data)):
  bandpass(vis=data[vis],caltable=data[vis]+'.bandpass.bcal',
        field=bpcal[vis],spw='',combine='scan',refant=refant,
        solint='inf',solnorm=T,minblperant=4,fillgaps=gaps,
        gaintable=[data[vis]+'.bpphase.gcal'])

for vis in range(len(data)):
  bandpass(vis=data[vis],caltable=data[vis]+'.bandpass_bpoly.bcal',
        field=bpcal[vis],spw='',combine='scan',refant=refant,
        solint='inf',solnorm=T,minblperant=4,fillgaps=gaps,
        bandtype='BPOLY',degamp=7,degphase=7,
        gaintable=[data[vis]+'.bpphase.gcal'])

To plot all our tables we will use both our AnalysisUtils and plotcal. The aU.plotbandpass command will create plot files for each combination of dataset and antenna, for both amplitude and phase. Inspect all these plots to make sure that the bandpass solutions look good. In Figures 6 and 7 we show an output sample for this command for the two bandpass calibrators used -- note how noisy the weaker bandpass calibrator (J1924) is without BPOLY.

# Set some plotting things
SPW=['0','1','2','3']
numants=16 # max for any of the input datasets
os.system('rm -rf cal_plots')
os.system('mkdir cal_plots')

os.system('rm -rf cal_plots/*bandpass_bpoly.bcal.*png')
for vis in range(len(data)):
  aU.plotbandpass(caltable=data[vis]+'.bandpass.bcal',
                caltable2=data[vis]+'.bandpass_bpoly.bcal',
                field=bpcal[vis],xaxis='freq',yaxis='both',
                figfile='cal_plots/'+data[vis]+'.bandpass_bpoly.png',
                interactive=False,subplot=42)

The next step in the calibration is to calculate amplitude and phase gains vs time for our calibrators. The ideal case here is to have a solution per integration of the data, but in this case we will need to use a solution interval of 30 seconds to avoid having many failed solutions, especially in the weak phase calibrator. First, we will calculate gains for phase and later using that information, we will solve for amplitude and phase. In all the next three executions we are using the bandpass calibration table, as it provides the gains for phase and amplitude vs frequency.

#Using 30s (5 integrations) per solution to avoid many failed solution of the weak calibrator.
for vis in range(len(data)):
  gaincal(vis=data[vis],caltable=data[vis]+'.intphase.gcal',
        field=calfields[vis],spw=calchan,refant=refant,
        calmode='p',solint='30s',minsnr=2.0,minblperant=4,
        gaintable=[data[vis]+'.bandpass_bpoly.bcal'])

To plot our *.intphase.gcal tables we use the next command, which will create files containing the phase gains vs time for all the antennas and for all the datasets. In Figure 7 we show an example of such plots. Again, you will need to check all plots to make sure the solutions are good.

Fig. 8. Phase gain solutions vs time for 5 antennas for all calibrators in the case of X90c.
# Plotting phase gains vs time
os.system('rm -rf cal_plots/*intphase*png')
for vis in data:
  for spw in SPW:
    for antenna in range(0,numants,5):    
      plotcal(caltable=vis+'.intphase.gcal',
        xaxis='time',yaxis='phase',antenna='%d~%d'%(antenna,antenna+4),
        iteration='antenna',subplot=511,poln='',spw=spw,
        showgui=F,
        figfile='cal_plots/'+vis+'.intphase.gcal.spw%s.ant%d_%d.png'%(spw,antenna,antenna+4),
        fontsize=8.0,plotrange=[0,0,-180,180])

Next we will use gaincal to solve for gain phases for all the calibrators, but this time we will get one single solution per scan.

# Gaincal execution as before, but to get a single solution per scan
for vis in range(len(data)):
  gaincal(vis=data[vis],caltable=data[vis]+'.scanphase.gcal',
        field=calfields[vis],spw=calchan,refant=refant,
        calmode='p',solint='inf',minsnr=2.0,minblperant=4,
        gaintable=[data[vis]+'.bandpass_bpoly.bcal'])

The next command will produce many plots, and like the previous one, you will get one for each dataset and for each antenna. See example in Figure 8.

Fig. 9. Phase gain solutions vs time for all the calibrators in the case of X90c showing the same 5 antennas from previous figure. Note that there is only one solution per scan for each source.
# Phase vs time plotting for our calibrators
os.system('rm -rf cal_plots/*scanphase*png')
for vis in data:
  for spw in SPW:
    for antenna in range(0,numants,5):
      plotcal(caltable=vis+'.scanphase.gcal',
        xaxis='time',yaxis='phase',antenna='%d~%d'%(antenna,antenna+4),
        iteration='antenna',subplot=511,poln='',spw=spw,
        showgui=F,
        figfile='cal_plots/'+vis+'.scanphase.gcal.spw%s.ant%d_%d.png'%(spw,antenna,antenna+4),
        fontsize=8.0,plotrange=[0,0,-180,180])

Finally, in the next gaincal we will solve for amplitude and phase. We will use the gain phase calibration table produced before. We will get one solution per scan for all our calibrators.

# Gaincal solving for amplitude and phase vs time.
for vis in range(len(data)):
  gaincal(vis=data[vis],caltable=data[vis]+'.amp.gcal',
        field=calfields[vis],spw=calchan,refant=refant,
        calmode='ap',solint='inf',minsnr=2.0,minblperant=4,
        gaintable=[data[vis]+'.bandpass_bpoly.bcal',data[vis]+'.intphase.gcal'])

We now can check the resulting plots with the next plotcal executions. In Figure 9 we show an example of these plots.

Fig. 10. Gain amplitudes for all the calibrators, as function of time. A plot for X90c is shown with 5 antennas in it. As in the previous plot, there is only one solution per scan for each source.
# Plotting gain amplitudes as function of time.
os.system('rm -rf cal_plots/*amp*png')
for vis in data:  
  for spw in SPW:
    for antenna in range(0,numants,5):
      plotcal(caltable=vis+'.amp.gcal',
        xaxis='time',yaxis='amp',antenna='%d~%d'%(antenna,antenna+4),
        iteration='antenna',subplot=511,poln='',spw=spw,
        figfile='cal_plots/'+vis+'.amp.gcal.spw%s.ant%d_%d.png'%(spw,antenna,antenna+4),
        fontsize=8.0)

We can plot the same tables but in a different way that will allow us to look for higher abnormal gains in the solutions. With the next command we will get only four plots for each spw. Each of those plots shows all the gains for all the antennas for all the sources (See Figure 10 for an example). Below in the box, there are some comments to focus you on some data, so you can double check. Before continuing, make sure you check all the calibration tables.

Fig. 11. Gain amplitudes for all the calibrators, as function of time (X90c), all antennas for each spw. Check for higher gains that deviate from the average.
# Look for low or high gains compared to other data
os.system('rm -rf *amp.png')
for vis in data:
  plotcal(caltable=vis+'.amp.gcal',
        xaxis='time',yaxis='amp',antenna='',field='',
        iteration='spw',subplot=411,poln='',spw='',
        showgui=False,figfile=vis+'.amp.png')
# X90c good
# X575 one antenna spw=0
# Xb50 end times bad all spws (low el), spw=0 more; low DV02 on Juno
# X39b low DV02 on Juno

Set Absolute Flux Scale

Now that the gain calibration is done, we need to set the flux for our calibrators. For this we will use Juno, our primary flux calibrator. We will do this by using fluxscale. We then will transfer the flux information from our phase calibrator to our science target. For these datasets, we note the following for this step: The derived flux densities for X39b, X90c, and Xb50 are quite reasonable. The results for X575 are a little high, probably because Juno was at low elevation for these observations.

# Setting fluxes

for vis in range(len(data)):
  fluxscale(vis=data[vis],caltable=data[vis]+'.amp.gcal',
            fluxtable=data[vis]+'.flux.cal',reference=fluxcal)

Below we copy the fluxscale numbers from the logger window:

#:calibrater::open     Opening MS: uid___A002_X3d4118_X39b.antwvrtsys.ms.fixed for calibration.
#:Calibrater:: Initializing nominal selection to the whole MS.
#:calibrater::fluxscale        Beginning fluxscale--(MSSelection version)-------
#:::    Assuming all non-reference fields are transfer fields.
#:::    Found reference field(s): Juno
#:::    Found transfer field(s):  1924-292 nrao530 ph 1625-254
#:::    Flux density for 1924-292 in SpW=0 is: 2.29951 +/- 0.0615168 (SNR = 37.3802, N = 24)
#:::    Flux density for 1924-292 in SpW=1 is: 2.3884 +/- 0.0706186 (SNR = 33.8212, N = 24)
#:::    Flux density for 1924-292 in SpW=2 is: 2.30261 +/- 0.0561279 (SNR = 41.0244, N = 24)
#:::    Flux density for 1924-292 in SpW=3 is: 2.27951 +/- 0.0682752 (SNR = 33.3871, N = 20)
#:::    Flux density for nrao530 ph in SpW=0 is: 0.627546 +/- 0.0366209 (SNR = 17.1363, N = 24)
#:::    Flux density for nrao530 ph in SpW=1 is: 0.64101 +/- 0.0406892 (SNR = 15.7538, N = 24)
#:::    Flux density for nrao530 ph in SpW=2 is: 0.615775 +/- 0.0324313 (SNR = 18.9871, N = 24)
#:::    Flux density for nrao530 ph in SpW=3 is: 0.607186 +/- 0.0407003 (SNR = 14.9184, N = 20)
#:::    Flux density for 1625-254 in SpW=0 is: 0.405301 +/- 0.0270957 (SNR = 14.9581, N = 24)
#:::    Flux density for 1625-254 in SpW=1 is: 0.414875 +/- 0.0274183 (SNR = 15.1313, N = 24)
#:::    Flux density for 1625-254 in SpW=2 is: 0.411922 +/- 0.0212975 (SNR = 19.3413, N = 24)
#:::    Flux density for 1625-254 in SpW=3 is: 0.406868 +/- 0.0302855 (SNR = 13.4344, N = 20)
#:Calibrater::fluxscale        Storing result in uid___A002_X3d4118_X39b.antwvrtsys.ms.fixed.flux.cal
#:::   Writing solutions to table: uid___A002_X3d4118_X39b.antwvrtsys.ms.fixed.flux.cal

#:calibrater::open     Opening MS: uid___A002_X3d55cb_X575.antwvrtsys.ms.fixed for calibration.
#:Calibrater:: Initializing nominal selection to the whole MS.
#:calibrater::fluxscale        Beginning fluxscale--(MSSelection version)-------
#:::    Assuming all non-reference fields are transfer fields.
#:::    Found reference field(s): Juno
#:::    Found transfer field(s):  3c279 nrao530 ph 1625-254
#:::    Flux density for 3c279 in SpW=0 is: 10.0511 +/- 0.170046 (SNR = 59.1085, N = 22)
#:::    Flux density for 3c279 in SpW=1 is: 11.0567 +/- 0.204718 (SNR = 54.0092, N = 22)
#:::    Flux density for 3c279 in SpW=2 is: 10.8403 +/- 0.190514 (SNR = 56.9002, N = 22)
#:::    Flux density for 3c279 in SpW=3 is: 10.9256 +/- 0.213188 (SNR = 51.2484, N = 22)
#:::    Flux density for nrao530 ph in SpW=0 is: 0.895265 +/- 0.0469496 (SNR = 19.0686, N = 22)
#:::    Flux density for nrao530 ph in SpW=1 is: 0.962985 +/- 0.057688 (SNR = 16.693, N = 22)
#:::    Flux density for nrao530 ph in SpW=2 is: 0.95216 +/- 0.063454 (SNR = 15.0055, N = 22)
#:::    Flux density for nrao530 ph in SpW=3 is: 0.953193 +/- 0.0667595 (SNR = 14.278, N = 22)
#:::    Flux density for 1625-254 in SpW=0 is: 0.633491 +/- 0.054662 (SNR = 11.5892, N = 22)
#:::    Flux density for 1625-254 in SpW=1 is: 0.635812 +/- 0.0564569 (SNR = 11.2619, N = 22)
#:::    Flux density for 1625-254 in SpW=2 is: 0.62294 +/- 0.0523949 (SNR = 11.8893, N = 22)
#:::    Flux density for 1625-254 in SpW=3 is: 0.640096 +/- 0.0646176 (SNR = 9.90592, N = 22)
#:Calibrater::fluxscale        Storing result in uid___A002_X3d55cb_X575.antwvrtsys.ms.fixed.flux.cal
#:::   Writing solutions to table: uid___A002_X3d55cb_X575.antwvrtsys.ms.fixed.flux.cal

#:calibrater::open     Opening MS: uid___A002_X3d55cb_X90c.antwvrtsys.ms.fixed for calibration.
#:Calibrater:: Initializing nominal selection to the whole MS.
#:calibrater::fluxscale        Beginning fluxscale--(MSSelection version)-------
#:::    Assuming all non-reference fields are transfer fields.
#:::    Found reference field(s): Juno
#:::    Found transfer field(s):  1924-292 nrao530 ph 1625-254
#:::    Flux density for 1924-292 in SpW=0 is: 2.26482 +/- 0.0804822 (SNR = 28.1406, N = 22)
#:::    Flux density for 1924-292 in SpW=1 is: 2.30424 +/- 0.0688451 (SNR = 33.4698, N = 22)
#:::    Flux density for 1924-292 in SpW=2 is: 2.31636 +/- 0.063015 (SNR = 36.7589, N = 22)
#:::    Flux density for 1924-292 in SpW=3 is: 2.30091 +/- 0.0628643 (SNR = 36.6011, N = 22)
#:::    Flux density for nrao530 ph in SpW=0 is: 0.655383 +/- 0.0392507 (SNR = 16.6974, N = 22)
#:::    Flux density for nrao530 ph in SpW=1 is: 0.663148 +/- 0.0358838 (SNR = 18.4804, N = 22)
#:::    Flux density for nrao530 ph in SpW=2 is: 0.660628 +/- 0.0318434 (SNR = 20.7462, N = 22)
#:::    Flux density for nrao530 ph in SpW=3 is: 0.656224 +/- 0.0347546 (SNR = 18.8817, N = 22)
#:::    Flux density for 1625-254 in SpW=0 is: 0.423542 +/- 0.0340076 (SNR = 12.4543, N = 22)
#:::    Flux density for 1625-254 in SpW=1 is: 0.413005 +/- 0.0278793 (SNR = 14.8141, N = 22)
#:::    Flux density for 1625-254 in SpW=2 is: 0.438477 +/- 0.030868 (SNR = 14.2049, N = 22)
#:::    Flux density for 1625-254 in SpW=3 is: 0.434589 +/- 0.028559 (SNR = 15.2173, N = 22)
#:Calibrater::fluxscale        Storing result in uid___A002_X3d55cb_X90c.antwvrtsys.ms.fixed.flux.cal
#:::   Writing solutions to table: uid___A002_X3d55cb_X90c.antwvrtsys.ms.fixed.flux.cal

#:calibrater::open     Opening MS: uid___A002_X3d55cb_Xb50.antwvrtsys.ms.fixed for calibration.
#:Calibrater:: Initializing nominal selection to the whole MS.
#:calibrater::fluxscale        Beginning fluxscale--(MSSelection version)-------
#:::    Assuming all non-reference fields are transfer fields.
#:::    Found reference field(s): Juno
#:::    Found transfer field(s):  1924-292 nrao530 ph 1625-254
#:::    Flux density for 1924-292 in SpW=0 is: 2.39168 +/- 0.0655984 (SNR = 36.4594, N = 22)
#:::    Flux density for 1924-292 in SpW=1 is: 2.33502 +/- 0.0478833 (SNR = 48.7649, N = 22)
#:::    Flux density for 1924-292 in SpW=2 is: 2.26212 +/- 0.053858 (SNR = 42.0016, N = 22)
#:::    Flux density for 1924-292 in SpW=3 is: 2.28193 +/- 0.0551113 (SNR = 41.4059, N = 22)
#:::    Flux density for nrao530 ph in SpW=0 is: 0.686101 +/- 0.0515145 (SNR = 13.3186, N = 22)
#:::    Flux density for nrao530 ph in SpW=1 is: 0.696972 +/- 0.0392324 (SNR = 17.7652, N = 22)
#:::    Flux density for nrao530 ph in SpW=2 is: 0.683672 +/- 0.0422611 (SNR = 16.1773, N = 22)
#:::    Flux density for nrao530 ph in SpW=3 is: 0.69109 +/- 0.0405247 (SNR = 17.0535, N = 22)
#:::    Flux density for 1625-254 in SpW=0 is: 0.504527 +/- 0.0496797 (SNR = 10.1556, N = 22)
#:::    Flux density for 1625-254 in SpW=1 is: 0.482058 +/- 0.0293865 (SNR = 16.4041, N = 22)
#:::    Flux density for 1625-254 in SpW=2 is: 0.445422 +/- 0.0365173 (SNR = 12.1976, N = 22)
#:::    Flux density for 1625-254 in SpW=3 is: 0.458493 +/- 0.0345567 (SNR = 13.2678, N = 22)
#:Calibrater::fluxscale        Storing result in uid___A002_X3d55cb_Xb50.antwvrtsys.ms.fixed.flux.cal
#:::   Writing solutions to table: uid___A002_X3d55cb_Xb50.antwvrtsys.ms.fixed.flux.cal

When derived fluxes are too high and nothing else appears wrong with the data, the cause (especially at Band 9) is likely to be decorrelation. So we will favor the average of the lower values to explicitly set the flux densities based on the fluxscale results. We set the flux of the two bandpass calibrators for the convenience of having a fully calibrated dataset. The only ones that really matter are 1625 as the gain calibrator and the check source nrao530. Based on this, we proceed as follows to make the changes that are needed:

# Fixing the fluxes
flux1924=[2.3,0,0,0]
flux1625=[0.43,0,0,0]
fluxnrao530=[0.66,0,0,0]
flux3c279=[10.5,0,0,0]

datawith1924=['uid___A002_X3d4118_X39b.antwvrtsys.ms.fixed',
      'uid___A002_X3d55cb_X90c.antwvrtsys.ms.fixed',
      'uid___A002_X3d55cb_Xb50.antwvrtsys.ms.fixed']
for vis in datawith1924:
  setjy(vis=vis,field='1924-292',fluxdensity=flux1924,usescratch=F)
  setjy(vis=vis,field='1625-254',fluxdensity=flux1625,usescratch=F)
  setjy(vis=vis,field='nrao530*',fluxdensity=fluxnrao530,usescratch=F)

datawith3c279=['uid___A002_X3d55cb_X575.antwvrtsys.ms.fixed']
for vis in datawith3c279:
  setjy(vis=vis,field='3c279',fluxdensity=flux3c279,usescratch=F)
  setjy(vis=vis,field='1625-254',fluxdensity=flux1625,usescratch=F)
  setjy(vis=vis,field='nrao530*',fluxdensity=fluxnrao530,usescratch=F)

We now need to re-run the amplitude calibration step.

#This new amplitude calibration will be used in the applycal.

data=['uid___A002_X3d4118_X39b.antwvrtsys.ms.fixed',
      'uid___A002_X3d55cb_X575.antwvrtsys.ms.fixed',
      'uid___A002_X3d55cb_X90c.antwvrtsys.ms.fixed',
      'uid___A002_X3d55cb_Xb50.antwvrtsys.ms.fixed']
# Match up intents with source names
pcal='1625-254'
fluxcal='Juno' 
science='IRAS16293*'
check='nrao530*'
bpcal=['1924-292','3c279','1924-292','1924-292']
calfields=['1924-292,Juno,1625-254,nrao530*',
           '3c279,Juno,1625-254,nrao530*',
           '1924-292,Juno,1625-254,nrao530*',
           '1924-292,Juno,1625-254,nrao530*']

calchan='0~3:20~3820'
refant='DV14'

for vis in range(len(data)):
  gaincal(vis=data[vis],caltable=data[vis]+'.amp.final.gcal',
        field=calfields[vis],spw=calchan,refant=refant,
        calmode='ap',solint='inf',minsnr=2.0,minblperant=4,
        gaintable=[data[vis]+'.bandpass_bpoly.bcal',data[vis]+'.intphase.gcal'])

Application of calibration tables

Now that we have all the calibration tables with the results, both for gains and flux, we need to apply all the tables to the data. First, we save the state of the flags.

data=['uid___A002_X3d4118_X39b.antwvrtsys.ms.fixed',
      'uid___A002_X3d55cb_X575.antwvrtsys.ms.fixed',
      'uid___A002_X3d55cb_X90c.antwvrtsys.ms.fixed',
      'uid___A002_X3d55cb_Xb50.antwvrtsys.ms.fixed']
# Match up intents with source names
pcal='1625-254'
fluxcal='Juno' 
science='IRAS16293*'
check='nrao530*'
bpcal=['1924-292','3c279','1924-292','1924-292']

for vis in range(len(data)):
  flagmanager(vis=data[vis],mode='save',versionname='beforeapplycal')

for vis in range(len(data)):
  applycal(vis=data[vis],field=bpcal[vis],
        gaintable=[data[vis]+'.bandpass_bpoly.bcal',data[vis]+'.intphase.gcal',
                   data[vis]+'.amp.final.gcal'],
        interp=['nearest','nearest','nearest'],
        gainfield=[bpcal[vis],bpcal[vis],bpcal[vis]],flagbackup=F,calwt=F)

for vis in range(len(data)):
  applycal(vis=data[vis],field=pcal,
         gaintable=[data[vis]+'.bandpass_bpoly.bcal',data[vis]+'.intphase.gcal',
                    data[vis]+'.amp.final.gcal'],
        interp=['nearest','nearest','nearest'],
         gainfield=[bpcal[vis],pcal,pcal],flagbackup=F,calwt=F)

for vis in range(len(data)):
  applycal(vis=data[vis],field=fluxcal,
         gaintable=[data[vis]+'.bandpass_bpoly.bcal',data[vis]+'.intphase.gcal',
                    data[vis]+'.amp.final.gcal'],
        interp=['nearest','nearest','nearest'],
         gainfield=[bpcal[vis],fluxcal,fluxcal],flagbackup=F,calwt=F)

for vis in range(len(data)):
  applycal(vis=data[vis],field=science,
        interp=['nearest','linear','linear'],
        gaintable=[data[vis]+'.bandpass_bpoly.bcal',data[vis]+'.scanphase.gcal',
                   data[vis]+'.amp.final.gcal'],
        gainfield=[bpcal[vis],pcal,pcal],flagbackup=F,calwt=F)

for vis in range(len(data)):
  applycal(vis=data[vis],field=check,
        interp=['nearest','linear','linear'],
        gaintable=[data[vis]+'.bandpass_bpoly.bcal',data[vis]+'.scanphase.gcal',
                   data[vis]+'.amp.final.gcal'],
        gainfield=[bpcal[vis],pcal,pcal],flagbackup=F,calwt=F)


Plot corrected data

Fig. 11. This plot shows amplitude (flux) versus time for the dataset X90c in spw 1.
Fig. 12. Phase vs time for all the sources in X90c (spw 1).
Fig. 13. Amplitude (flux) vs frequency for 1625-254 in X39b.
Fig. 14. Phase vs frequency for 3c279 in X575.
Fig. 15. Amplitude vs frequency for the science target, colored by fields for X39b and spw 1, showing CO (6-5).
Fig. 16. Amplitude vs frequency as for Figure 15, but for spw 3.

The next commands will help you visualize the result of the application of the calibration tables to the data. You can check if the amplitudes and phases vs time and frequency look reasonable for all the sources, in particular for the science target fields.

The next command will produce four plots, as the one we showed in Figure 4, but the amplitude in this new plots corresponds to flux because we now have calibrated data. See Figure 11 for an example of it. It is important to check that all the sources have similar amplitude (flux) in the different spws and datasets.

data=['uid___A002_X3d4118_X39b.antwvrtsys.ms.fixed',
      'uid___A002_X3d55cb_X575.antwvrtsys.ms.fixed',
      'uid___A002_X3d55cb_X90c.antwvrtsys.ms.fixed',
      'uid___A002_X3d55cb_Xb50.antwvrtsys.ms.fixed']
# Match up intents with source names
pcal='1625-254'
fluxcal='Juno' 
science='IRAS16293*'
check='nrao530*'
bpcal=['1924-292','3c279','1924-292','1924-292']
calfields=['1924-292,Juno,1625-254,nrao530*',
           '3c279,Juno,1625-254,nrao530*',
           '1924-292,Juno,1625-254,nrao530*',
           '1924-292,Juno,1625-254,nrao530*']
# Set some plotting things
SPW=['0','1','2','3']

os.system('rm -rf aftercal_plots')
os.system('mkdir aftercal_plots')

os.system('rm -rf aftercal_plots/*cal.time*.png')
for vis in range(len(data)):
  for spw in SPW:
    plotms(vis=data[vis],spw=spw,xaxis='time',yaxis='amp',field='',avgchannel='3840',
       coloraxis='field',ydatacolumn='corrected',
       plotfile='aftercal_plots/'+data[vis]+'.cal.time.spw%s.png'%(spw))

It is also important to check the phases vs time for all the sources. The next command will get you the corresponding plots, four for each dataset. You can see that while the bandpass and amplitude calibrator have very concentrated phases around 0 degrees, the phase calibrator and the science target do not. In Figure 12 we show an example of these plots for spw 3 in the dataset X90c.

os.system('rm -rf aftercal_plots/*cal.time.phase*.png')
for vis in range(len(data)):
  for spw in SPW:
    plotms(vis=data[vis],spw=spw,xaxis='time',yaxis='phase',field='',avgchannel='3840',
       coloraxis='field',ydatacolumn='corrected',
       plotfile='aftercal_plots/'+data[vis]+'.cal.time.phase.spw%s.png'%(spw))

We now check the amplitude of the sources vs frequency. This is important since we expect that all the spw have very similar behavior. You will have one plot for each field for each dataset. In Figure 13 we show the case for 1625-254, our phase calibrator, for X90c.

os.system('rm -rf aftercal_plots/*cal.freq.amp*.png')
for vis in range(len(data)):
  for field in calfields[vis].split(','):
    plotms(vis=data[vis],field='%s'%field,xaxis='freq', yaxis='amp',
       spw='',avgtime='1e8',avgscan=T,
       coloraxis='spw',xselfscale=T,ydatacolumn='corrected',
       plotfile='aftercal_plots/'+data[vis]+'.cal.freq.amp.'+field+'.png')

The next command will produce similar plots but this time of phase vs frequency (see Figure 14 for an example). You will notice that only strong sources, like 3c279, will show clearly phases concentrated around 0 degrees.

os.system('rm -rf aftercal_plots/*cal.freq.phase*.png')
for vis in range(len(data)):
  for field in calfields[vis].split(','):
    plotms(vis=data[vis],field='%s'%field,xaxis='freq', yaxis='phase',
       spw='',avgtime='1e8',avgscan=T,
       coloraxis='spw',xselfscale=T,ydatacolumn='corrected',
       plotfile='aftercal_plots/'+data[vis]+'.cal.freq.phase.'+field+'.png')

Finally, for our science target, we plot amplitude (flux) vs frequency for all the spectral windows. See an example of this is Figure 15 and 16. If you look carefully at these plots, you will notice that the datasets at low elevation show much more line emission - this is due to the shorter projected baselines present for these datasets (X39b and Xb50). In contrast, you start to see some weak absorption for the dataset at high elevation: X90c.

os.system('rm -rf aftercal_plots/*science.freq.amp*.png')
for vis in range(len(data)):
  for spw in SPW:    
    plotms(vis=data[vis],field=science,xaxis='freq', yaxis='amp',
       spw=spw,avgtime='1e8',avgscan=T,
       coloraxis='field',xselfscale=T,ydatacolumn='corrected',
       plotfile='aftercal_plots/'+data[vis]+'.science.freq.amp'+spw+'.png',
       title=data[vis].split('_')[-1].split('.')[0]+'.IRAS16293_spw'+spw)

Next, we put a list of additional plotms commands that do not produce .png files, but that you can explore and save a copy of the output file if you are interested.

# Additional manual plots
vis=data[0]
plotms(vis=vis,spw='',xaxis='time',yaxis='amp',field='',avgchannel='3840',
       coloraxis='field',ydatacolumn='corrected',iteraxis='spw')

vis=data[1]
plotms(vis=vis,field=bpcal[1],xaxis='freq', yaxis='amp',
       spw='2',avgtime='1e8',avgscan=T,iteraxis='antenna',
       coloraxis='field',xselfscale=T,ydatacolumn='corrected')

vis=data[1]
plotms(vis=vis,field=calfields[1],xaxis='freq', yaxis='amp',
       spw='',avgtime='1e8',avgscan=T,
       coloraxis='field',xselfscale=T,ydatacolumn='corrected')

Split and concatenate the calibrated data

You are ready now to extract the final calibrated data for the science target. Of course, you can do something similar for the calibrators if you are interested.

# Splitting final calibrated datasets
data=['uid___A002_X3d4118_X39b.antwvrtsys.ms.fixed',
      'uid___A002_X3d55cb_X575.antwvrtsys.ms.fixed',
      'uid___A002_X3d55cb_X90c.antwvrtsys.ms.fixed',
      'uid___A002_X3d55cb_Xb50.antwvrtsys.ms.fixed']

for vis in data:
  split(vis=vis,outputvis='%s.cal.IRAS16293.ms.fixed'%(vis.split('.')[0]),
     datacolumn='corrected',field='IRAS16293*',keepflags=False)

There are now four datasets fully calibrated. We finally will merge these datasets into a single one, in order to proceed with the analysis and imaging in both continuum and spectral line.

# Concatenating the final split files
concatdata=['uid___A002_X3d4118_X39b.cal.IRAS16293.ms.fixed',
     'uid___A002_X3d55cb_X575.cal.IRAS16293.ms.fixed',
     'uid___A002_X3d55cb_X90c.cal.IRAS16293.ms.fixed',
     'uid___A002_X3d55cb_Xb50.cal.IRAS16293.ms.fixed']

concat(vis=concatdata,concatvis='IRAS16293_Band9.fixed.ms')

To speed up imaging we apply 60 second time averaging to the concatenated dataset.

# 60s time averaging
split(vis='IRAS16293_Band9.fixed.ms', datacolumn='data', timebin='60s', 
       outputvis='IRAS16293_Band9.fixed.rebin.ms')

As a final step, we zero the rows of the pointing table because it is quite large and is not currently needed by the imaging software for mosaics -indeed its presence will cause an error during imaging if you skip this step.

# Remove rows of pointing table
tb.open('IRAS16293_Band9.fixed.rebin.ms/POINTING', nomodify = False)
a = tb.rownumbers()
tb.removerows(a)
tb.close()

Now you have completed the calibration and have everything you need to carry out the imaging stage. Follow IRAS16293 Band9 - Imaging to go to the imaging section of this casaguide.

Last checked on CASA Version 3.4.0.