Difference between revisions of "EVLA high frequency spectral line tutorial - IRC+10216 - calibration"

From CASA Guides
Jump to: navigation, search
(Set Up the Model for the Flux Calibrator)
(Analyze the Image Cubes)
 
(9 intermediate revisions by the same user not shown)
Line 103: Line 103:
  
 
Create a plot of antenna positions using {{plotants}}.
 
Create a plot of antenna positions using {{plotants}}.
 
[[Image:elevationvstime.png|thumb|Elevation as a function of time (after selecting colorize by field).]]
 
  
 
<source lang="python">
 
<source lang="python">
Line 178: Line 176:
 
The absolute fluxes for the frequencies have now been determined and one can proceed to the bandpass and complex gain calibrations.
 
The absolute fluxes for the frequencies have now been determined and one can proceed to the bandpass and complex gain calibrations.
  
==Bandpass==
+
==Bandpass Calibration==
 
 
Before determining the bandpass solution, we need to inspect phase and amplitude
 
variations with time and frequency on the bandpass calibrator to
 
decide how best to proceed. We limit the number of antennas to make
 
the plot easier to see. We chose ea02 as it seems like a good
 
candidate for the reference antenna.
 
 
 
<source lang="python">
 
# In CASA
 
plotms(vis=vis,field='5',
 
      xaxis='channel',yaxis='phase',correlation='RR',
 
      avgtime='1e8',spw='0:4~60',antenna='ea02&ea23')
 
</source>
 
[[Image:Nobandpass_phase.png|thumb|Phase as a function of channel for ea02 (after Custom and upping "Style" to 3.)]]
 
The phase variation is modest ~10 degrees.  Now expand to all baselines that include ea02, then hit the "Plot" button.
 
 
 
<source lang="python">
 
# In CASA
 
plotms(vis=vis,field='5',
 
      xaxis='channel',yaxis='phase',correlation='RR',
 
      avgtime='1e8',spw='0:4~60',antenna='ea02', coloraxis='antenna2')
 
</source>
 
[[Image:Nobandpass_phasetime.png|thumb|Phase as a function of time for all baselines with antenna ea02 (after Custom and upping "Style" to 3.)]]
 
From this
 
you can see that the phase variation across the bandpass is
 
modest. Next check LL, and spw=1, both correlations. Also check
 
other antennas if you like.
 
 
 
Now look at the phase as a function of time.
 
 
 
<source lang="python">
 
# In CASA
 
plotms(vis=vis,field='5',
 
      xaxis='time',yaxis='phase',correlation='RR',
 
      avgchannel='64',spw='0:4~60',antenna='ea02&ea23')
 
</source>
 
 
 
  
Expand to all antennas with ea02
+
Correct the phase variations with time before solving for the bandpass to
 
 
<source lang="python">
 
# In CASA
 
plotms(vis=vis,field='5',
 
      xaxis='time',yaxis='phase',correlation='RR',
 
      avgchannel='64',spw='0:4~60',antenna='ea02', coloraxis='antenna2')
 
</source>
 
 
 
You may want to select "Custom" under "unflagged points symbol" and then change Style from "2" to "3" under the "Display" tab, "Unflagged Points Symbol" .
 
 
 
You can see that the phase variations are smooth, but do vary
 
significantly over the 5 minutes of observation -- in most cases by
 
a few 10s of degrees. Zoom in to see this better if you want.
 
 
 
The conclusion from this investigation is that you need to correct
 
the phase variations with time before solving for the bandpass to
 
 
prevent decorrelation of the vector averaged bandpass
 
prevent decorrelation of the vector averaged bandpass
solution. Since the phase variation as a function of channel is
+
solution.  
modest, you can average over several channels to increase the signal
 
to noise of the phase vs. time solution. If the phase variation as a
 
function of channel is larger you may need to use only a few
 
channels to prevent introducing delay-based closure errors as can happen from averaging over
 
non-bandpass corrected channels with large phase variations.
 
  
 +
[[Image:Prebp_phasecal2.png|thumb|Phase only calibration before bandpass. The 4 lines are both polarizations in both spw, unfortunately two of them get the same color green at the moment.]]
  
Since the bandpass calibrator is quite strong we do the phase-only
 
solution on the integration time of 20 seconds (solint='int').
 
 
[[Image:Prebp_phasecal2.png|thumb|Phase only calibration before bandpass. The 4 lines are both polarizations in both spw, unfortunately two of them get the same color green at the moment.]]
 
 
<source lang="python">
 
<source lang="python">
 
# In CASA
 
# In CASA
gaincal(vis=vis,caltable='bpphase.gcal',
+
gaincal(vis=vis, caltable='bpphase.gcal', field='5', spw='0~1:20~40',
        field='5',spw='0~1:20~40',
+
         refant='ea02', gaintype='G', calmode='p',solint='int',
         refant='ea02',calmode='p',solint='int',minsnr=2.0,
 
 
         gaincurve=T)
 
         gaincurve=T)
 
</source>
 
</source>
  
Plot the solutions
+
Plot the solutions using {{plotcal}}.
  
 
<source lang="python">
 
<source lang="python">
Line 262: Line 198:
 
         iteration='antenna',subplot=331,plotrange=[0,0,-180,180])
 
         iteration='antenna',subplot=331,plotrange=[0,0,-180,180])
 
</source>
 
</source>
 
These solutions will appear in the CASA plotter gui. If you closed it after plotting the antennas above, it should reopen. If it is still open from before, the new plots should just appear. After you are done looking at the first set of plots, push the "Next" button on the GUI to see the next set of antennas.
 
  
 
Next we can apply this phase solution on the fly while determining
 
Next we can apply this phase solution on the fly while determining
Line 270: Line 204:
 
<source lang="python">
 
<source lang="python">
 
# In CASA
 
# In CASA
bandpass(vis=vis,caltable='bandpass.bcal',field='5',
+
bandpass(vis=vis, caltable='bandpass.bcal', field='5', refant='ea02',
        refant='ea02',solint='inf',solnorm=T,
+
        solint='inf', solnorm=T, gaintable='bpphase.gcal',
        gaintable='bpphase.gcal',
+
        gaincurve=T)
        gaincurve=T)
 
 
</source>
 
</source>
 
'''A few words about solint and combine:'''
 
 
The use of  solint='inf' in {{bandpass}} will derive one bandpass
 
solution for the whole J1229+0203 scan. Note that if there had been two observations of the bandpass calibrator (for example), this command would have combined the data from both scans to form one bandpass solution, because the default of the combine parameter '''for {{bandpass}}''' is combine='scan'. To solve for one bandpass for each bandpass calibrator scan you would also need to include combine='''' '''' in the bandpass call. In all calibration tasks, regardless of solint, scan boundaries are only crossed when combine='scan'. Likewise, field (spw) boundaries are only crossed if combine='field' (combine='spw'), the latter two are not generally good ideas for bandpass solutions. 
 
  
 
Plot the solutions, amplitude and phase:
 
Plot the solutions, amplitude and phase:
Line 298: Line 226:
 
</source>
 
</source>
  
Note that phases on ea12 look noiser than on other antennas. This
+
==Remaining Calibration==
jitter could indicate bad pointing; note that ea12 had just come back in
 
the array.
 
 
 
This step isn't necessary from a calibration perspective, but if you
 
want to go ahead and check the bandpass calibration on the bandpass
 
calibrator you can run {{applycal}} here. In future we hope to plot
 
corrected data on-the-fly without this {{applycal}} step. Later applycals
 
will overwrite this one, so no need to worry.
 
 
 
[[Image:Applybandpass_phase.png|thumb|Phase as a function of channel, plotting the corrected data (after Custom and upping "Style" to 3.)]]
 
 
 
<source lang="python">
 
applycal(vis=vis,field='5',
 
        gaintable='bandpass.bcal',
 
        gainfield=['','5'],
 
        gaincurve=T,calwt=F)
 
</source>
 
 
 
Similar to {{flagdata}}, {{applycal}} works like a matrix. The first entries in the lists are to be used together, so are the second entries etc. All will be applied to the 'field' selection. In the above example, the 'bandpass.cal' that was obtained for field '5' (the bandpass observation) is also applied to field '2'. Again, lists within the lists are fine. 
 
 
 
<source lang="python">
 
plotms(vis=vis,field='5',
 
      xaxis='channel',yaxis='phase',ydatacolumn='corrected',
 
      correlation='RR',
 
      avgtime='1e8',spw='0:4~60',antenna='ea02', coloraxis='antenna2')
 
</source>
 
 
 
<source lang="python">
 
plotms(vis=vis,field='5',
 
      xaxis='channel',yaxis='amp',ydatacolumn='corrected',
 
      correlation='RR',
 
      avgtime='1e8',spw='0:4~60',antenna='ea02', coloraxis='antenna2')
 
</source>
 
 
 
Note that the phase and amplitude as a function of channel are very flat now.
 
 
 
==Gain Calibration==
 
 
 
Now that we have a bandpass solution to apply we can solve for the antenna-based phase and amplitude gain calibration. Since the phase changes on a much shorter timescale than the amplitude, we will solve for them separately. In particular, if the phase changes significantly over a scan time, the amplitude would be decorrelated, if the un-corrected phase were averaged over this timescale. Note that we re-solve for the gain solutions of the bandpass calibrator, so we can derive new solutions that are corrected for the bandpass shape. Since the bandpass calibrator will not be used again, this is not strictly necessary, but is useful to check its calibrated flux density for example. We use a minimum signal-to-noise of 2 here as it seems to be a good compromise for using good data without rejecting too many solutions (minsnr=2).
 
 
 
<source lang="python">
 
# In CASA
 
gaincal(vis=vis,caltable='intphase.gcal',
 
        field='2,5,7',spw='0~1:4~60',
 
        refant='ea02',calmode='p',solint='int',minsnr=2.0,
 
        gaintable='bandpass.bcal',
 
        gaincurve=T)
 
</source>
 
[[Image:allcal_phaseint2.png|thumb|Plot of phase solutions on an integration time.]]
 
 
 
Here solint='int' coupled with calmode='p' will derive a single phase solution for each 10 second integration. Note that the bandpass table is applied on-the-fly before solving for the phase solutions, however the bandpass is NOT applied to the data permanently until applycal is run later on.
 
 
 
Note that quite a few solutions are rejected due to SNR<2 (printed to terminal). For the most part it
 
is only one or two solutions out of >30 so this isn't too worrying. Take note if you see large numbers of rejected solutions per integration. This is likely an indication that solint is too short for the S/N of the data.
 
 
 
Now look at the phase solution, and note the obvious scatter within a scan time.
 
 
 
<source lang="python">
 
# In CASA
 
plotcal(caltable='intphase.gcal',xaxis='time',yaxis='phase',
 
        iteration='antenna',subplot=331,plotrange=[0,0,-180,180])
 
</source>
 
 
 
Although solint='int' (i.e. the integration time of 20 seconds) is the best choice to apply before for solving for the amplitude solutions, it is not a good idea to use this to apply to the target. This is because the phase-scatter within a scan can dominate the interpolation between calibrator scans. Instead, we also solve for the phase on the scan time, solint='inf' (but combine='''' '''', since we want one solution per scan) for application to the target later on. '''Unlike the bandpass task,''' for gaincal, the default of the combine parameter is combine='''' ''''.
 
[[Image:allcal_phaseinf2.png|thumb|Plot of phase solutions on a scan time.]]
 
<source lang="python">
 
# In CASA
 
gaincal(vis=vis,caltable='scanphase.gcal',
 
        field='2,5,7',spw='0~1:4~60',
 
        refant='ea02',calmode='p',solint='inf',minsnr=2.0,
 
        gaintable='bandpass.bcal',
 
        gaincurve=T)
 
</source>
 
 
 
<source lang="python">
 
# In CASA
 
plotcal(caltable='scanphase.gcal',xaxis='time',yaxis='phase',
 
        iteration='antenna',subplot=331,plotrange=[0,0,-180,180])
 
</source>
 
 
 
Note that there are no failed solutions here because of the added S/N afforded by the longer solint.
 
Alternatively, instead of making a separate phase solution for application to the target, one can also run {{smoothcal}} to smooth the solutions derived on the integration time.
 
 
 
Next we apply the bandpass and solint='int' phase-only calibration solutions on-the-fly to derive amplitude solutions.
 
Here the use of solint='inf', but combine='''' '''' will result in one solution per scan interval.
 
 
 
<source lang="python">
 
# In CASA
 
gaincal(vis=vis,caltable='amp.gcal',
 
        field='2,5,7',spw='0~1:4~60',
 
        refant='ea02',calmode='ap',solint='inf',minsnr=2.0,
 
        gaintable=['bandpass.bcal','intphase.gcal'],
 
        gaincurve=T)
 
</source>
 
[[Image:allcal_ampphase.png|thumb|Plot of residual phase solutions on a scan time]]
 
 
 
Now let's look at the resulting phase solutions.  Since we have taken out the phase as best we can by applying the solint='int' phase-only solution, this plot will give a good idea of the residual phase error. If you see scatter of more than a few degrees here, you should consider going back and looking for more data to flag, particularly bad timeranges etc.
 
 
 
<source lang="python">
 
# In CASA     
 
plotcal(caltable='amp.gcal',xaxis='time',yaxis='phase',
 
        iteration='antenna',subplot=331)
 
</source>
 
 
 
Indeed, both antenna ea12 (all times) and ea23 (first 1/3 of observation) show particularly large residual phase noise.
 
[[Image:allcal_amp.png|thumb|Plot of amplitude solutions on a scan time]]
 
<source lang="python">
 
# In CASA
 
plotcal(caltable='amp.gcal',xaxis='time',yaxis='amp',
 
        iteration='antenna',subplot=331)
 
</source>
 
 
 
Note that the amplitude solutions for ea12 are very low; this is another indication that this antenna is dubious.
 
 
 
Next we use the flux calibrator (whose flux density was set in {{setjy}} above) to derive the flux of the other calibrators. Note that the flux table REPLACES the amp.gcal in terms of future application of the calibration to the data, i.e. the flux table contains both the amp.gcal and flux scaling. Unlike the gain calibration steps, this is not an incremental table.
 
 
 
<source lang="python">
 
# In CASA
 
fluxscale(vis=vis,caltable='amp.gcal',
 
          fluxtable='flux.cal',reference='7')
 
</source>
 
 
 
[[Image:allcal_flux.png|thumb|Plot of flux corrected amplitude solutions.]]
 
It is a good idea to note down for your records the derived flux densities:
 
 
 
<pre style="background-color: #fffacd;">
 
Flux density for J0954+1743 in SpW=0 is: 0.237135 +/- 0.000858511 (SNR = 276.216, nAnt= 19)
 
Flux density for J0954+1743 in SpW=1 is: 0.247597 +/- 0.00063498 (SNR = 389.928, nAnt= 19)
 
Flux density for J1229+0203 in SpW=0 is: 26.508 +/- 0 (SNR = inf, nAnt= 19)
 
Flux density for J1229+0203 in SpW=1 is: 26.2517 +/- 0 (SNR = inf, nAnt= 19)
 
 
 
</pre>
 
 
 
Obviously, the signal-to-noise for J1229+0203 can't be infinity! This is just an indication that there is only one scan for this source, and we derived a scan based amplitude solution, so there is no variation to calculate.
 
 
 
Next, check that the flux.cal table looks as expected.
 
 
 
<source lang="python">
 
# In CASA
 
plotcal(caltable='flux.cal',xaxis='time',yaxis='amp',
 
        iteration='antenna',subplot=331)
 
</source>
 
 
 
==Applycal and Inspect==
 
 
 
Now we apply the calibration to each source, according to which tables are appropriate, and which source should be used to do that particular calibration. For the calibrators, all bandpass solutions come from the bandpass calibrator (id=5), and the phase and amplitude calibration comes from their own solutions.
 
 
 
'''Note:''' In applycal we set calwt=F. It is very important to turn off this parameter which determines if the weights are calibrated along with the data. Data from antennas with better receiver performance and/or longer integration times should have higher weights, and it can be advantageous to factor this information into the calibration. During the VLA era, meaningful weights were available for each visibility. However, EVLA is not yet recording the information necessary to calculate meaningful weights. Since these data weights are used at the imaging stage you can get strange results from having calwt=T when the input weights are themselves not meaningful, especially for self-calibration on resolved sources (your flux calibrator and target for example). In a few months EVLA data will again have meaningful weights and the default calwt=T will likely again be the best option.
 
 
 
<source lang="python">
 
# In CASA
 
# for the gain/phase calibrator (field '2'):
 
applycal(vis=vis,field='2',
 
        gaintable=['bandpass.bcal','intphase.gcal','flux.cal'],
 
        gainfield=['','5','2','2'],
 
        gaincurve=T,calwt=F)
 
</source>
 
 
 
<source lang="python">
 
# In CASA
 
# for the bandpass calibrator (field '5'):
 
applycal(vis=vis,field='5',
 
        gaintable=['bandpass.bcal','intphase.gcal','flux.cal'],
 
        gainfield=['','5','5','5'],
 
        gaincurve=T,calwt=F)
 
</source>
 
 
 
<source lang="python">
 
# In CASA
 
# for the flux calibrator (field '7'):
 
applycal(vis=vis,field='7',
 
        gaintable=['bandpass.bcal','intphase.gcal','flux.cal'],
 
        gainfield=['','5','7','7'],
 
        gaincurve=T,calwt=F)
 
</source>
 
 
 
For the target we apply the bandpass from id=5, and the calibration from the gain calibrator (id=2):
 
  
<source lang="python">
+
In this Community Day Event guide, we will skip over the remaining calibration steps.  However, you can refer to [[EVLA high frequency Spectral Line tutorial - IRC+10216 part1]] for full details. To summarize the missing steps, you bootstrap the flux densities of the secondary calibrators by
# In CASA
 
# for the target source IRC10216 (field '3'):
 
applycal(vis=vis,field='3',
 
        gaintable=['bandpass.bcal','scanphase.gcal','flux.cal'],
 
        gainfield=['','5','2','2'],
 
        gaincurve=T,calwt=F)
 
</source>
 
  
Now inspect the corrected data:
+
# doing a phase only calibration on all calibrators using {{gaincal}},
[[Image:applycal_inspect.png|thumb|Plot of calibrated amplitudes over time.]]
+
# doing amplitude only calibration on all calibrators while applying the phase-only solutions, and
<source lang="python">
+
# deriving the flux density of the secondary calibrators while applying the previously obtained solutions.
# In CASA
 
plotms(vis=vis,field='5',ydatacolumn='corrected',
 
      xaxis='time',yaxis='amp',correlation='RR,LL',
 
      avgchannel='64',spw='0:4~60',antenna='', coloraxis='antenna1')
 
</source>
 
 
 
This plot shows some data deviating from the average amplitudes. Use methods described above to
 
mark a region for a small number of deviant data points, and click "Locate". You will find that ea12 is responsible.
 
  
<source lang="python">
+
Then, you can calibrate the target source using the phase and amplitude solutions you have obtained.  Apply the target source calibration solutions using {{applycal}}.  Use {{plotms}} to examine the calibrated data.  If more flagging is required, redo all calibration steps.  When the data look good, {{split}} the target source into a separate measurement set.  Subtract the continuum flux using {{uvcontsub}}.  Make Doppler corrections using {{cvel}}, or let {{clean}} do the Doppler corrections on the fly.
# In CASA
 
plotms(vis=vis,field='2',ydatacolumn='corrected',
 
      xaxis='time',yaxis='amp',correlation='RR,LL',
 
      avgchannel='64',spw='0:4~60',antenna='', coloraxis='antenna2')
 
</source>
 
  
Here we see some problems, with high points. Mark some regions
+
This tutorial picks up where [[ EVLA high frequency spectral line tutorial - IRC+10216 - calibration]] leaves off.  
and locate in {{plotms}} to find out which antennas and in which spws. Pay special
 
attention to antennas that have been called out already as showing some dubious behavior.
 
  
What you find is that ea07 which we flagged spw=1 above, is also bad for the same timerange in spw=0. This was not obvious in the raw data, because spw=0 was adjusted in the on-line system by a gain attenuator, while spw=1 wasn't. So a lack of power on this antenna can look like very low (and obvious) amplitudes in spw=1 but not for spw=0. Looking carefully you'll see that ea07 is actually pretty noisy throughout.
+
==Imaging==
[[Image:ea12.png|thumb|Plot of antenna ea12 by itself]]
 
[[Image:ea23.png|thumb|Plot of antenna ea23 by itself]]
 
  
From the locate we also find that ea12 and ea23 show some high points; to see this, replot baselines using each of them alone:
+
[[Image:irc10216_uvspec.png|thumb|UV-plot of the spectral line signal in both spw for IRC+10216.]]
  
<source lang="python">
+
The continuum-subtracted spectral line data is contained in <tt>IRC10216_spls.ms</tt>.  Use {{plotms}} to plot the lines.
plotms(vis=vis,field='2',ydatacolumn='corrected',
 
      xaxis='time',yaxis='amp',correlation='RR,LL',
 
      avgchannel='64',spw='0:4~60',antenna='ea12', coloraxis='antenna2')
 
</source>
 
  
 
<source lang="python">
 
<source lang="python">
plotms(vis=vis,field='2',ydatacolumn='corrected',
+
plotms(vis='IRC10216_spls.ms', xaxis='channel', yaxis='amp',
      xaxis='time',yaxis='amp',correlation='RR,LL',
+
       avgtime='1e8', avgscan=T, coloraxis='spw')
       avgchannel='64',spw='0:4~60',antenna='ea23', coloraxis='antenna2')
 
 
</source>
 
</source>
  
it may be a a good idea to flag ea12 completely - it's just a bit noisy all around and ea23 is pretty noisy during the first scans between initial and second pointing. Recall that these are antennas we became suspicious of while inspecting the calibration solutions.
+
Now it is time to image the visibility data using {{clean}}. For illustration, we will clean channel 22 of the SiS line.
 
 
[[Image:target_uvdist.png|thumb|IRC+10216 as a function of uv-distance.]]
 
Now lets see how the target looks. Because the target has resolved structure, its best to look at it as
 
a function of uvdistance. We'll go ahead and exclude the three antennas we already know have problems.
 
  
 
<source lang="python">
 
<source lang="python">
 
# In CASA
 
# In CASA
plotms(vis=vis,field='3',ydatacolumn='corrected',
+
os.system('rm -rf ch22.*') # remove previously generated image, if it exists
      xaxis='uvdist',yaxis='amp',correlation='RR,LL',
+
clean(vis='IRC10216_spls.ms', imagename='ch22', spw='1:22~22',
      avgchannel='64',spw='0:4~60',antenna='!ea07;!ea12;!ea23', coloraxis='antenna2')
+
      mode='channel', nchan=1, start='', width=1, niter=100000,
 +
      gain=0.1, threshold='3.0mJy', psfmode='clark', imagermode='csclean',
 +
      interactive=T, npercycle=100, imsize=300, cell=['0.4arcsec', '0.4arcsec'],
 +
      stokes='I', weighting='briggs', robust=0.5)  
 
</source>
 
</source>
  
you can see that the spikes
+
After running the above command, make a region in the CASA viewer.  Double-click inside the region and clean by clicking the green circular arrow. After each cycle, click the green circular arrow again if the flux inside the region is brighter than the flux peaks outside the region.  When you are finished cleaning, click the red 'STOP' button.  
are caused by a single antenna. Use, zoom, mark, and locate to see which one.
 
Also look at spw=1.
 
  
Turns out to be ea28; to confirm, replot with antenna=!ea28, and
+
Open the resulting image using the {{viewer}}.
  
 
<source lang="python">
 
<source lang="python">
 
# In CASA
 
# In CASA
plotms(vis=vis,field='3',ydatacolumn='corrected',
+
viewer("ch22.image")
      xaxis='uvdist',yaxis='amp',correlation='RR,LL',
 
      avgchannel='64',spw='0:4~60',antenna='!ea07;!ea12;!ea23;!ea28', coloraxis='antenna2')
 
 
</source>
 
</source>
  
To see if it's restricted to a certain time
+
In the viewer, make a region and click inside it to get statistics about the portion of the image in the region.
  
<source lang="python">
+
==Analyze the Image Cubes==
# In CASA
 
plotms(vis=vis,field='3',ydatacolumn='corrected',
 
      xaxis='time',yaxis='amp',correlation='RR,LL',
 
      avgchannel='64',spw='0:4~60',antenna='ea28', coloraxis='antenna1')
 
</source>
 
  
 
+
Open the HC3N image cube in the {{viewer}}.
Baselines with ea28 clearly show issues until about two-thirds of the way through the observation.
 
Plot another distant antenna to compare. We will go ahead and flag it all, since its hanging far out on the north
 
arm by itself.
 
 
 
The additional data we've identified as bad need to be flagged, and then all the calibration steps will need to be run
 
again.
 
  
 
<source lang="python">
 
<source lang="python">
 
# In CASA
 
# In CASA
flagdata(vis=vis,
+
viewer("IRC10216_HC3N.image")
        field=['',''],
 
        spw=['',''],
 
        antenna=['ea07,ea12,ea28','ea07,ea23'],
 
        timerange=['','03:21:40~04:10:00'])
 
 
</source>
 
</source>
  
==Redo Calibration after more Flagging==
+
Click the 'Play' button to view all channels in the cube.  Select '''Tool --> Spectral profile...''' from the pull down menu.  Create a region on the image to view the spectrum within the region.
  
After flagging, you'll need to repeat the calibration steps above.  Here, we append  _redo to the table names to distinguish them from the first round, in case we want to compare with previous versions.  
+
Open the SiS image cube in the viewer.
  
 
<source lang="python">
 
<source lang="python">
 
# In CASA
 
# In CASA
gaincal(vis=vis,caltable='bpphase_redo.gcal',
+
viewer("IRC10216_SiS.image")
        field='5',spw='0~1:20~40',
 
        refant='ea02',calmode='p',solint='int',minsnr=2.0,
 
        gaincurve=T)
 
 
</source>
 
</source>
  
<source lang="python">
+
Determine what channels in the cube have emission. Then make moment maps using {{immoments}}.
# In CASA
 
bandpass(vis=vis,caltable='bandpass_redo.bcal',
 
        field='5',
 
        refant='ea02',solint='inf',solnorm=T,
 
        gaintable=['bpphase_redo.gcal'],
 
        gaincurve=T)
 
</source>
 
  
 
<source lang="python">
 
<source lang="python">
 
# In CASA
 
# In CASA
gaincal(vis=vis,caltable='intphase_redo.gcal',
+
os.system('rm -rf IRC10216_Sis.mom0') # remove previously generated map, if it exists
        field='2,5,7',spw='0~1:4~60',
+
immoments(imagename="IRC10216_SiS.image", moments=[0], axis="spectral",
        refant='ea02',calmode='p',solint='int',minsnr=2.0,
+
          chans="12~40", outfile="IRC10216_Sis.mom0")  
        gaintable=['bandpass_redo.bcal'],
 
        gaincurve=T)
 
 
</source>
 
</source>
  
<source lang="python">
+
Open the moment map in the viewer.
# In CASA
 
gaincal(vis=vis,caltable='scanphase_redo.gcal',
 
        field='2,5,7',spw='0~1:4~60',
 
        refant='ea02',calmode='p',solint='inf',minsnr=2.0,
 
        gaintable=['bandpass_redo.bcal'],
 
        gaincurve=T)
 
</source>
 
  
 
<source lang="python">
 
<source lang="python">
 
# In CASA
 
# In CASA
gaincal(vis=vis,caltable='amp_redo.gcal',
+
viewer("IRC10216_Sis.mom0")
        field='2,5,7',spw='0~1:4~60',
 
        refant='ea02',calmode='ap',solint='inf',minsnr=2.0,
 
        gaintable=['bandpass_redo.bcal','intphase_redo.gcal'],
 
        gaincurve=T)
 
 
</source>
 
</source>
 
<source lang="python">
 
# In CASA     
 
fluxscale(vis=vis,caltable='amp_redo.gcal',
 
          fluxtable='flux_redo.cal',reference='7')
 
</source>
 
 
<pre style="background-color: #fffacd;">
 
Flux density for J0954+1743 in SpW=0 is: 0.235345 +/- 0.000879422
 
    (SNR = 267.613, nAnt= 16)
 
Flux density for J0954+1743 in SpW=1 is: 0.241996 +/- 0.000930228
 
    (SNR = 260.147, nAnt= 16)
 
Flux density for J1229+0203 in SpW=0 is: 25.2479 +/- 0
 
    (SNR = inf, nAnt= 16)
 
Flux density for J1229+0203 in SpW=1 is: 24.9907 +/- 0
 
    (SNR = inf, nAnt= 16)
 
  
Flux density for J0954+1743 in SpW=0 is: 0.247052 +/- 0.000946345 (SNR = 261.059, nAnt= 16)
+
Overlay contours by selecting '''Data --> Open''' from the pull down menu, selecting the moment-0 map, and clicking 'contour map'.
Flux density for J0954+1743 in SpW=1 is: 0.254038 +/- 0.00097531 (SNR = 260.47, nAnt= 16)
 
Flux density for J1229+0203 in SpW=0 is: 26.5079 +/- 0 (SNR = inf, nAnt= 16)
 
Flux density for J1229+0203 in SpW=1 is: 26.2335 +/- 0 (SNR = inf, nAnt= 16)
 
</pre>
 
 
 
Feel free to pause here and remake the calibration solution plots from above, just be sure to put in the revised table names.
 
 
 
==Redo Applycal and Inspect==
 
 
 
Now, apply all the new calibrations, which will overwrite the old ones.  These commands are identical to those above, with the exception of the _redo part of each calibration filename.
 
 
 
<source lang="python">
 
# In CASA
 
applycal(vis=vis,field='2',
 
        gaintable=['bandpass_redo.bcal','intphase_redo.gcal','flux_redo.cal'],
 
        gainfield=['','5','2','2'],
 
        gaincurve=T,calwt=F)
 
</source>
 
 
 
<source lang="python">
 
# In CASA
 
applycal(vis=vis,field='5',
 
        gaintable=['bandpass_redo.bcal','intphase_redo.gcal','flux_redo.cal'],
 
        gainfield=['','5','5','5'],
 
        gaincurve=T,calwt=F)
 
</source>
 
 
 
<source lang="python">
 
# In CASA
 
applycal(vis=vis,field='7',
 
        gaintable=['bandpass_redo.bcal','intphase_redo.gcal','flux_redo.cal'],
 
        gainfield=['','5','7','7'],
 
        gaincurve=T,calwt=F)
 
</source>
 
 
 
[[Image:gaincal_corrflag.png|thumb|Gain calibrator after further flagging and recalibration]]
 
[[Image:target_corrflag.png|thumb|IRC+10216 after further flagging and recalibration (after selecting colorize by spw).]]
 
 
 
<source lang="python">
 
# In CASA
 
applycal(vis=vis,field='3',
 
        gaintable=['bandpass_redo.bcal','scanphase_redo.gcal','flux_redo.cal'],
 
        gainfield=['','5','2','2'],
 
        gaincurve=T,calwt=F)
 
</source>
 
 
 
Now you can inspect the calibrated data again. Except for random scatter things look pretty good.
 
 
 
<source lang="python">
 
# In CASA
 
plotms(vis=vis,field='2',ydatacolumn='corrected',
 
      xaxis='time',yaxis='amp',correlation='RR,LL',
 
      avgchannel='64',spw='0:4~60',antenna='', coloraxis='antenna2')
 
</source>
 
 
 
You can use the Mark and Locate buttons to assess that the remaining scatter seems random, i.e. no particular antenna or time range appears to be responsible.
 
 
 
<source lang="python">
 
# In CASA
 
plotms(vis=vis,field='3',ydatacolumn='corrected',
 
      xaxis='uvdist',yaxis='amp',correlation='RR,LL',
 
      avgchannel='64',spw='0~1:4~60',antenna='', coloraxis='spw')
 
</source>
 
 
 
==Split==
 
 
 
Now we split the data into individual files. This is not strictly necessary, as you can select the appropriate fields in later clean stages, but it is safer in case for example you get confused with later processing and want to fall back to this point (this is especially a good idea if you plan to do continuum subtraction or self calibration later on). It also makes smaller individual files in case you want to copy to another machine or colleague.
 
 
 
Here, we split off the data for the phase calibrator and the target:
 
 
 
<source lang="python">
 
# In CASA
 
split(vis=vis,outputvis='J0954',
 
      field='2')
 
</source>
 
 
 
<source lang="python">
 
# In CASA
 
split(vis=vis,outputvis='IRC10216',
 
      field='3')
 
</source>
 
 
 
To reinitialize the scratch columns for use by later tasks, we need to run clearcal for both new datasets
 
 
 
<source lang="python">
 
# In CASA
 
clearcal(vis='J0954')
 
</source>
 
 
 
<source lang="python">
 
# In CASA
 
clearcal(vis='IRC10216')
 
</source>
 
 
 
This concludes the calibration phase of the data reductions. The tutorial continues with continuum subtraction, imaging, and image analysis in
 
[[EVLA high frequency spectral line tutorial - IRC+10216 - imaging]].
 
  
 
[[Main Page | &#8629; '''CASAguides''']]
 
[[Main Page | &#8629; '''CASAguides''']]
 
--[[User:Cbrogan|Crystal Brogan]]   
 
--additions:[[User:Jott| Juergen Ott]]
 
  
 
{{Checked 3.3.0}}
 
{{Checked 3.3.0}}

Latest revision as of 16:57, 6 January 2012

This page is under development

This CASA Guide provides a simplified approach to the editing, calibration, and imaging of a galactic EVLA observation for use in NRAO Community Day Events. For a complete description of the data reduction process see EVLA high frequency Spectral Line tutorial - IRC+10216 part1 and EVLA high frequency Spectral Line tutorial - IRC+10216 part2.

Overview

VLT V-band image of IRC+10216 showing dust rings out to a radius of 90" by Leão et al. (2006, A&A, 455, 187).

This tutorial describes the data reduction for two spectral lines, HC3N and SiS, observed toward the AGB star IRC+10216 by the EVLA in D-configuration.

Getting the data

The data for this tutorial can be obtained by anonymous FTP from ftp://ftp.aoc.nrao.edu/staff/gvanmoor/community_day/. Download all 4 TAR files.

For example,

# In UNIX
wget 'ftp://ftp.aoc.nrao.edu/staff/gvanmoor/community_day/*'

Initial Inspection and Flagging

GEt a summary listing of the data set using listobs.

# In CASA
vis = 'day2_TDEM0003_20s_full'
listobs(vis=vis, verbose=True)

Below we have cut and pasted the most relevant output from the logger.

Fields: 4
  ID   Code Name                RA              Decl          Epoch   SrcId nVis
  2    D    J0954+1743          09:54:56.82363 +17.43.31.2224 J2000   2     32726
  3    NONE IRC+10216           09:47:57.38200 +13.16.40.6600 J2000   3     99540
  5    F    J1229+0203          12:29:06.69973 +02.03.08.5982 J2000   5     5436
  7    E    J1331+3030          13:31:08.28798 +30.30.32.9589 J2000   7     2736
   (nVis = Total number of time/baseline visibilities per field)
Spectral Windows:  (2 unique spectral windows and 1 unique polarization setups)
  SpwID  #Chans Frame Ch1(MHz)    ChanWid(kHz)  TotBW(kHz)  Corrs
  0          64 TOPO  36387.2295  125           8000        RR  RL  LR  LL
  1          64 TOPO  36304.542   125           8000        RR  RL  LR  LL
Sources: 10
  ID   Name                SpwId RestFreq(MHz)  SysVel(km/s)
  0    J1008+0730          0     0.03639232     -0.026
  0    J1008+0730          1     0.03639232     -0.026
  2    J0954+1743          0     0.03639232     -0.026
  2    J0954+1743          1     0.03639232     -0.026
  3    IRC+10216           0     0.03639232     -0.026
  3    IRC+10216           1     0.03639232     -0.026
  5    J1229+0203          0     0.03639232     -0.026
  5    J1229+0203          1     0.03639232     -0.026
  7    J1331+3030          0     0.03639232     -0.026
  7    J1331+3030          1     0.03639232     -0.026
Antennas: 19:
  ID   Name  Station   Diam.    Long.         Lat.
  0    ea01  W09       25.0 m   -107.37.25.2  +33.53.51.0
  1    ea02  E02       25.0 m   -107.37.04.4  +33.54.01.1
  2    ea03  E09       25.0 m   -107.36.45.1  +33.53.53.6
  3    ea04  W01       25.0 m   -107.37.05.9  +33.54.00.5
  4    ea05  W08       25.0 m   -107.37.21.6  +33.53.53.0
  5    ea07  N06       25.0 m   -107.37.06.9  +33.54.10.3
  6    ea08  N01       25.0 m   -107.37.06.0  +33.54.01.8
  7    ea09  E06       25.0 m   -107.36.55.6  +33.53.57.7
  8    ea12  E08       25.0 m   -107.36.48.9  +33.53.55.1
  9    ea15  W06       25.0 m   -107.37.15.6  +33.53.56.4
  10   ea19  W04       25.0 m   -107.37.10.8  +33.53.59.1
  11   ea20  N05       25.0 m   -107.37.06.7  +33.54.08.0
  12   ea21  E01       25.0 m   -107.37.05.7  +33.53.59.2
  13   ea22  N04       25.0 m   -107.37.06.5  +33.54.06.1
  14   ea23  E07       25.0 m   -107.36.52.4  +33.53.56.5
  15   ea24  W05       25.0 m   -107.37.13.0  +33.53.57.8
  16   ea25  N02       25.0 m   -107.37.06.2  +33.54.03.5
  17   ea27  E03       25.0 m   -107.37.02.8  +33.54.00.5
  18   ea28  N08       25.0 m   -107.37.07.5  +33.54.15.8

We summarize the observing strategy in this table.

Gain calibrator J0954+1743 field id = 2
Bandpass calibrator J1229+0203 field id = 5
Flux calibrator J1331+3030 (3C286) field id = 7
Science target IRC+10216 field id = 3
Antenna locations from running plotants

Create a plot of antenna positions using plotants.

# In CASA
plotants(vis=vis)
Result of plotms
Zooming in and marking region (hatched box)

Next, let's look at all the source amplitudes as a function of time using plotms.

# In CASA
plotms(vis=vis, xaxis='time', yaxis='amp', correlation='RR,LL',
       avgchannel='64', spw='0:4~60', coloraxis='field')

Now zoom in on the region very near zero amplitude for sources J0954+1743 and IRC+10216. To zoom, select the Zoom tool in lower left corner of the plotms GUI, then you can left click to draw a box. Look for the low values (you may want to zoom a few times to really see the suspect points clearly). Now use the Mark Region and Locate buttons (located along the bottom of the GUI) to see which antenna is causing problems. The output will be shown in the logger. Since all the "located" baselines include ea12, this is the responsible antenna.

Now click the clear region button, and then go back to the zoom button to zoom in further to note exactly what the time range is: 03:41:00~04:10:00.

Check the other sideband by changing spw to 1:4~60. You will have to rezoom. If you have trouble, click on the Mark icon and then back to zoom. In spw=1, ea07 is bad from the beginning until after next pointing run: 03:21:40~04:10:00. To see this, compare the amplitudes when antenna is set to 'ea07' and when it is set to one of the other antennas, such as 'ea08'.

If you set antenna to 'ea12' and zoom in on this intial timerange, you can also see that ea12 is bad during the same time range as for spw 0. You can also see this by entering '!ea07' for antenna, which removes ea07 from the plot (in CASA selection, ! deselects).

We can flag the bad data using flagdata.

# In CASA
flagdata(vis=vis,
         field=['2,3','2,3'],
         spw=['','1'],
         antenna=['ea12','ea07'],
         timerange=['03:41:00~04:10:00','03:21:40~04:10:00'])

flagdata works by spanning up a matrix. The first entries in each list must be taken as one flagging command, as well as the second entries etc. Lists within lists are fine. In the above example, the first flagging command is issued for fields 2 and 3 for all spws and within the 03:41:00~04:10:00 timerange. A second command is again for the fields 2 and 3 but for spw 1 only and for the second timerange in the list '03:21:40~04:10:00'.

Note that because the chosen timerange is limited to fields 2 and 3, the field parameter is not really needed; however, flagdata will run fastest if you put as many constraints as possible.

Set Up the Model for the Flux Calibrator

We set the model for the flux calibrator using setjy. First, check the availability of calibration models.

# In CASA
setjy(vis=vis,listmodimages=T)

There is no Ka-band model of 3C286. We will use the K-band model instead.

# In CASA
setjy(vis=vis,field='7',spw='0~1',
      modimage='3C286_K.im')

setjy scales the total flux in the model image to that appropriate for your individual spectral window frequencies according to the calibrator's flux and reports this number to the logger.

The logger output for each spw is:
J1331+3030 (fld ind 7) spw 0  [I=1.7762, Q=0, U=0, V=0] Jy, (Perley-Butler 2010)
J1331+3030 (fld ind 7) spw 1  [I=1.7794, Q=0, U=0, V=0] Jy, (Perley-Butler 2010)

The absolute fluxes for the frequencies have now been determined and one can proceed to the bandpass and complex gain calibrations.

Bandpass Calibration

Correct the phase variations with time before solving for the bandpass to prevent decorrelation of the vector averaged bandpass solution.

Phase only calibration before bandpass. The 4 lines are both polarizations in both spw, unfortunately two of them get the same color green at the moment.
# In CASA
gaincal(vis=vis, caltable='bpphase.gcal', field='5', spw='0~1:20~40',
        refant='ea02', gaintype='G', calmode='p',solint='int',
        gaincurve=T)

Plot the solutions using plotcal.

# In CASA
plotcal(caltable='bpphase.gcal',xaxis='time',yaxis='phase',
        iteration='antenna',subplot=331,plotrange=[0,0,-180,180])

Next we can apply this phase solution on the fly while determining the bandpass solutions on the timescale of the bandpass calibrator scan (solint='inf').

# In CASA
bandpass(vis=vis, caltable='bandpass.bcal', field='5', refant='ea02',
         solint='inf', solnorm=T, gaintable='bpphase.gcal',
         gaincurve=T)

Plot the solutions, amplitude and phase:

Amplitude Bandpass solutions
Phase Bandpass solutions
# In CASA
plotcal(caltable='bandpass.bcal',xaxis='chan',yaxis='amp',
        iteration='antenna',subplot=331)


# In CASA
plotcal(caltable='bandpass.bcal',xaxis='chan',yaxis='phase',
        iteration='antenna',subplot=331)

Remaining Calibration

In this Community Day Event guide, we will skip over the remaining calibration steps. However, you can refer to EVLA high frequency Spectral Line tutorial - IRC+10216 part1 for full details. To summarize the missing steps, you bootstrap the flux densities of the secondary calibrators by

  1. doing a phase only calibration on all calibrators using gaincal,
  2. doing amplitude only calibration on all calibrators while applying the phase-only solutions, and
  3. deriving the flux density of the secondary calibrators while applying the previously obtained solutions.

Then, you can calibrate the target source using the phase and amplitude solutions you have obtained. Apply the target source calibration solutions using applycal. Use plotms to examine the calibrated data. If more flagging is required, redo all calibration steps. When the data look good, split the target source into a separate measurement set. Subtract the continuum flux using uvcontsub. Make Doppler corrections using cvel, or let clean do the Doppler corrections on the fly.

This tutorial picks up where EVLA high frequency spectral line tutorial - IRC+10216 - calibration leaves off.

Imaging

UV-plot of the spectral line signal in both spw for IRC+10216.

The continuum-subtracted spectral line data is contained in IRC10216_spls.ms. Use plotms to plot the lines.

plotms(vis='IRC10216_spls.ms', xaxis='channel', yaxis='amp',
       avgtime='1e8', avgscan=T, coloraxis='spw')

Now it is time to image the visibility data using clean. For illustration, we will clean channel 22 of the SiS line.

# In CASA
os.system('rm -rf ch22.*') # remove previously generated image, if it exists
clean(vis='IRC10216_spls.ms', imagename='ch22', spw='1:22~22',
      mode='channel', nchan=1, start='', width=1, niter=100000, 
      gain=0.1, threshold='3.0mJy', psfmode='clark', imagermode='csclean',
      interactive=T, npercycle=100, imsize=300, cell=['0.4arcsec', '0.4arcsec'],
      stokes='I', weighting='briggs', robust=0.5)

After running the above command, make a region in the CASA viewer. Double-click inside the region and clean by clicking the green circular arrow. After each cycle, click the green circular arrow again if the flux inside the region is brighter than the flux peaks outside the region. When you are finished cleaning, click the red 'STOP' button.

Open the resulting image using the viewer.

# In CASA
viewer("ch22.image")

In the viewer, make a region and click inside it to get statistics about the portion of the image in the region.

Analyze the Image Cubes

Open the HC3N image cube in the viewer.

# In CASA
viewer("IRC10216_HC3N.image")

Click the 'Play' button to view all channels in the cube. Select Tool --> Spectral profile... from the pull down menu. Create a region on the image to view the spectrum within the region.

Open the SiS image cube in the viewer.

# In CASA
viewer("IRC10216_SiS.image")

Determine what channels in the cube have emission. Then make moment maps using immoments.

# In CASA
os.system('rm -rf IRC10216_Sis.mom0') # remove previously generated map, if it exists
immoments(imagename="IRC10216_SiS.image", moments=[0], axis="spectral",
          chans="12~40", outfile="IRC10216_Sis.mom0")

Open the moment map in the viewer.

# In CASA
viewer("IRC10216_Sis.mom0")

Overlay contours by selecting Data --> Open from the pull down menu, selecting the moment-0 map, and clicking 'contour map'.

CASAguides

Last checked on CASA Version 3.3.0.