A script used for data reduction (2007-12-5)
#!/bin/csh -f
##############################################################
# Edmund Hodges-Kluck (UMd) & Jun-Hui Zhao (CfA) 12/5/07
# for the reduction of SMA data taken recently (2007 or later).
# The SMA data, used in this example, were taken in a spectral
# configuration with a hybrid resolution. The observations had
# multiple pointing for mosaic image of molecular lines. This
# E2E script consists of nine major steps with some detailed
# comments
# 1. READ THE RAW SMA DATA INTO MIRIAD
# 2. INSPECT THE RAW UV DATA
# 3. FRINGE AMPLITUDE CORRECTION
# 4. BANDPASS CALIBRATION
# 5. ANTENNA GAIN CALIBRATION
# 6. BOOTSTRAP THE FLUX DENSITY SCALE
# 7. CONTINUUM IMAGE
# 8. LINE/CONTINUUM SEPARATION
# 9. CONSTRUCT LINE IMAGE CUBE
##############################################################
### INTRODUCTION #############################################
##############################################################
# This is a full reduction in one line (HCO+ 3-2) looking at
# a data set with data in both HCO+ and HCN 3-2 in the lower
# sideband. To run the script, sequentially step through and
# uncomment the various 'continue' markers and run through
# until the script exits, then comment the marker out and
# uncomment the marker for the next step. Refer to the SMA
# MIRIAD manual for additional information. Note that we
# find a non-detection on the source, but get nice images of
# Sgr A*.
##############################################################
set source=S1
goto continue
### LOAD THE DATA ############################################
# In this section we take the data from its original MIR
# form and load it into MIRIAD.
##############################################################
#continue:
set darea=/2006/science/mir_data/071004_00:26:54/
set darea2=/2006/science/mir_data/071004_04:41:13/
echo "1. READ THE RAW SMA DATA INTO MIRIAD"
\rm -r 1E1740_1_rx1.lsb 1E1740_2_rx1.lsb
smalod in=$darea out=1E1740_1 rxif=1 \
sideband=0 readant=8 nscans=1
smalod in=$darea2 out=1E1740_2 rxif=1 \
sideband=0 readant=8 nscans=1
# We had 2 files for the observation, so load in the data
# twice and merge.
\rm -r 1E1740_rx1.lsb
uvaver vis=1E1740_1_rx1.lsb,1E1740_2_rx1.lsb out=1E1740_rx1.lsb
exit
### INSPECT RAW DATA ########################################
# In this section we look first at the sources in the
# observation (uvindex) and then the amplitude vs.
# time and the system temperature vs. time
#############################################################
#continue:
echo "2. INSPECT THE RAW UV DATA"
echo "List the observing scans in the uv data"
uvindex vis=1E1740_rx1.lsb
echo "Press to continue"
$<
smauvplt vis=1E1740_rx1.lsb \
axis=time,ampl device=/xs nxy=1,1 \
yrange=0,300
echo "Press to continue"
$<
smavarplt vis=1E1740_rx1.lsb device=/xs \
xaxis=time \
yaxis=systemp nxy=2,4
exit
### FRINGE AMPLITUDE CORRECTION #############################
# The changes in the system temperature over the duration
# of the observation essentially track changes in the
# atmosphere. The system temperature measurements are
# obtained by comparing the receiver output power from
# an ambient temperature load with the sky power, so we
# can calibrate the amplitude gains.
#############################################################
#continue:
echo "3. FRINGE AMPLITUDE CORRECTION"
\rm -r 1E1740_rx1.lsb.tsys
smafix vis=1E1740_rx1.lsb out=1E1740_rx1.lsb.tsys \
device=/xs xaxis=time \
yaxis=systemp nxy=2,4 \
options=tsyscorr
echo "Inspect vis data and flag bad data"
# We could filter out various times while we're doing the
# flagging here. We can also do this later when looking
# at the bandpass calibrator (as is the case in this script)
# Jupiter was just there for pointing and makes it harder to do
# interactive flagging.
uvflag vis=1E1740_rx1.lsb.tsys select='source(jupiter)' flagval=flag
# Antenna 4 was apparently not on, flag it just in case
uvflag vis=1E1740_rx1.lsb.tsys select='ant(4)' \
flagval=flag
# Interactively flag bad data
smablflag vis=1E1740_rx1.lsb.tsys axis=time,ampl device=/xs
exit
#continue:
# Look at the bandpass calibrator
echo "Inspect spectral data"
echo "for uranus"
smauvspec vis=1E1740_rx1.lsb.tsys \
select='source(uranus)' \
interval=1000 \
axis=freq,both device=/xs nxy=1,5 hann=5
# SgrA* too
#continue:
smauvspec vis=1E1740_rx1.lsb.tsys \
select='source(sgra*,Sg*)' \
interval=1000 \
axis=freq,both device=/xs nxy=1,5 hann=5
exit
### BANDPASS CALIBRATION ####################################
# Bandpass calibration gives the frequency solution, so
# we use something bright -- in the case of the SMA, a
# planet typically. In this case, Uranus.
#############################################################
#continue:
echo "4. BANDPASS CALIBRATION"
# smamfcal makes calibration corrections (antenna gains, delay
# terms and passband shapes) to a multi-frequency observation.
# The gains are calculated periodically using the interval
# specified. If phase stability is poor or you are using a
# planet, set weight=2 (see manual). weight=2 is default for
# the SMA.
smamfcal vis=1E1740_rx1.lsb.tsys select='source(uranus)' edge=2.2 \
weight=2 refant=3 interval=600
echo "Check up the bandpass solutions"
smagpplt vis=1E1740_rx1.lsb.tsys device=/xs \
yaxis=amp,phase \
options=bandpass polyfit=5 \
nxy=2,4
#exit
#continue:
echo 'apply the bpass solutions to the data'
\rm -r 1E1740_rx1.lsb.tsys.bp
uvaver vis=1E1740_rx1.lsb.tsys out=1E1740_rx1.lsb.tsys.bp
smauvspec vis=1E1740_rx1.lsb.tsys.bp \
select='source(sgr*,Sgr*)' \
interval=1000 \
axis=freq,both device=/xs nxy=2,2
exit
#continue:
# Interactively flag again to get anything still messing up
# the solution.
smablflag vis=1E1740_rx1.lsb.tsys.bp \
axis=time,ampl device=/xs
exit
#continue:
uvflag vis=1E1740_rx1.lsb.tsys.bp edge=15,15 \
flagval=flag select='window(1,2,3,4)'
# We have to vary the windows we're looking to clean up.
# For this observation,
# select='window(20,21,22,23,24)' selects the HCN (3-2)
# select='window(1,2,3,4)' selects the HCO+ (use this for now,
# it's brighter!)
smauvspec vis=1E1740_rx1.lsb.tsys.bp \
select='source(sgr*,Sgr*),ant(1)(6)' \
interval=1000 \
axis=freq,both device=/xs nxy=1,1
exit
### ANTENNA GAIN CALIBRATION ################################
# The gain calibration sources for this observation were
# 1733-130, SgrAstar. Here we apply the gain calibration to
# the data.
#############################################################
#continue:
echo "5. ANTENNA GAIN CALIBRATION"
# Filter out bad times here -- we decide if/what to filter
# based on the smauvplt below, so may have to run this step
# again with the line below uncommented and changed as
# appropriate.
# NOTE: May need to exclude 1.75h-3h based on the
# system temperature/phase problems.
#
# uvflag vis=1E1740_rx1.lsb.tsys.bp \
# select='time(1:45,3:00)' flagval=flag
#
# If you are unhappy with the result of the flagging, you can
# return to this step and run:
#
# uvflag vis=1E1740_rx1.lsb.tsys.bp \
# select='time(1:45,3:00)' flagval=unflag
#
# However, doing so may unflag other data you didn't want to unflag,
# so you might have to run smablflag again to interactively flag
# the data again. If you are worried, start over and don't flag
# based on time. You might also create a new file to do the time
# filtering on.
smablflag vis=1E1740_rx1.lsb.tsys.bp device=/xs axis=real,imag \
select='source(1733*,Sgr*)' options=nobase
echo "Press to continue."
$<
smauvplt vis=1E1740_rx1.lsb.tsys.bp \
select='source(1733*,Sgr*)' \
device=/xs options=nocal axis=time,amp nxy=1,3
echo "Split the gain calibrators from the multiple source file"
#continue:
# The number after these files will differ based on the dataset.
# Therefore, this step may need to be run twice.
\rm -r 1733-130.267745 SgrAstar.267745
uvsplit vis=1E1740_rx1.lsb.tsys.bp \
select='source(1733-130,SgrAstar)' options=nowindow
echo "rename the calibrator files"
\rm -r 1733-130.tsys.bp sgra.tsys.bp
mv 1733-130.267745 1733-130.tsys.bp
mv SgrAStar.267745 sgra.tsys.bp
exit
#continue:
echo " "
echo "Solve for antenna gains"
echo " "
# Can we use selfcal? Let's hope so. They're both strong point
# sources in this case.
# Note that here we do phase only so that we can solve for
# phase only.
selfcal vis=1733-130.tsys.bp interval=1 options=phas refant=3
selfcal vis=sgra.tsys.bp interval=1 options=phas refant=3
echo "Plot the gain solutions"
smagpplt vis=1733-130.tsys.bp device=/xs \
yaxis=phase log=0102+584gain.log \
options=gains polyfit=6 \
nxy=2,4
echo "Press to continue."
$<
smagpplt vis=sgra.tsys.bp device=/xs \
yaxis=phas \
options=gains polyfit=6 \
nxy=2,4
exit
#continue:
echo "Merge the two gain tables"
smagpplt vis=1733-130.tsys.bp,sgra.tsys.bp device=/xs \
yaxis=amp,pha \
options=gains,merge polyfit=8 \
nxy=2,4
#the merged gain table is place in the second file, namely,
#sgra.tsys.bp
smagpplt vis=sgra.tsys.bp device=/xs \
yaxis=phas \
options=gains polyfit=6 \
nxy=2,4
# Note that the default interpolation interval equals to
# the solution interval. One may make a longer interval
# in order to apply (interpolate) the selfcal solutions to
# the target sources.
# Modify the SgrA* header file so that the interval value
# is changed to 0.2day, then copy the gain table back to the
# main file with all sources.
puthd in=sgra.tsys.bp/interval value=0.2
gpcopy vis=sgra.tsys.bp out=1E1740_rx1.lsb.tsys.bp \
mode=copy options=nopass
# See what the results are.
smauvplt vis=1E1740_rx1.lsb.tsys.bp axis=time,amp device=/xs nxy=1,4
smauvspec vis=1E1740_rx1.lsb.tsys.bp axis=freq,both \
interval=1000000 options=nobase,avall \
select='source(sgra*)' device=/xs nxy=1,1 hann=15
echo 'apply the gain to the data'
\rm -r 1E1740_rx1.lsb.tsys.bp.gn
uvaver vis=1E1740_rx1.lsb.tsys.bp out=1E1740_rx1.lsb.tsys.bp.gn
exit
### BOOTSTRAP THE FLUX DENSITY SCALE ########################
# Find the scaling factor by looking at the bandpass
# calibrator, then apply the gain to the data
#############################################################
#continue:
echo "6. BOOTSTRAP THE FLUX DENSITY SCALE"
# This will calculate the scaling factor based on your
# bandpass calibrator (the SMA function looks specifically
# for a planet) and then applies it to the data. To
# calculate but not apply, use the option "noapply".
smaflux vis=1E1740_rx1.lsb.tsys.bp.gn \
select='source(uranus)' mirhome=$MIR
# Specific to this dataset:
# Found solar system objects: uranus
# Frequency(GHz) Tb(Kelvin)
# 267.797 93.61
# Scaling the data by 0.838
echo "Apply the scaling factor gains to the data"
\rm -r 1E1740_rx1.lsb.tsys.bp.gn.av
uvaver vis=1E1740_rx1.lsb.tsys.bp.gn out=1E1740_rx1.lsb.tsys.bp.gn.av
# Check for the spectrum of Sgr A*
smauvspec vis=1E1740_rx1.lsb.tsys.bp.gn.av \
select='source(sgra*),window(1,2,3,4)' \
interval=1000 hann=5 \
axis=freq,both device=/xs nxy=1,1 options=nobase,avall
exit
#############################################################
# At this point, the data has had all the calibrations applied
# to it. We now look at various sources within the data set.
#############################################################
# Change this to look at different sources, but remember to change
# manually in other places too as necessary
#continue:
set source=SgrAStar
### CONTINUUM IMAGE #########################################
# We make an image (via inversion) of the source in the
# continuum bands -- in this case we put lines in the 1,2,3,4
# and 20,21,22,23,24 channels and had the rest as wide channels
# with no lines. Thus, we look at the source in all of these
# wide channels.
#############################################################
echo "7. CONTINUUM IMAGE"
\rm -r $source.con.map $source.con.beam $source.con.icmp
invert vis=1E1740_rx1.lsb.tsys.bp.gn.av \
map=$source.con.map beam=$source.con.beam \
imsize=512 cell=0.6 sup=0 \
select='window(5,6,7,8,9,10,11,12,13,14,15,16,17,18,19),source(SgrA*)' \
line=chan,1,1,480
clean map=$source.con.map beam=$source.con.beam \
out=$source.con.icmp gain=0.08 \
cutoff=0 niters=3500 \
region=arcsec,'box(-25,-25,25,25)'
\rm -r $source.con.icln
restor model=$source.con.icmp beam=$source.con.beam \
map= $source.con.map \
out=$source.con.icln
smacgdisp in=$source.con.icln type=pixel \
region=arcsec,'boxes(-25,-25,25,25)' \
xybin=1,1 device=/xs nxy=1,1 \
options=full,beambr,wedge,trlab,3val \
labtyp=arcsec,arcsec range=0,0,lin,2 cols1=7 csize=1,1,1
imstat in=$source.con.icln region='box(20,20,80,80)'
exit
### LINE/CONTINUUM SEPARATION ################################
# In this section we subtract the continuum from the lines to
# enhance contrasts in what is left.
##############################################################
continue:
set source=S1
echo "8. LINE/CONTINUUM SEPARATION"
# In this case, our rest frequency was set to 267.55762 (HCO+ 3-2)
\rm -r $source.hco+.uv
# You must set source here -- $source won't work
uvaver vis=1E1740_rx1.lsb.tsys.bp.gn.av \
select='source(S1_CR),window(1,2,3,4)'\
out=$source.hco+.uv
# Plot of amplitude vs. LSR velocity and phase vs. velocity
smauvspec vis=$source.hco+.uv \
interval=10000 options=nobase,avall \
axis=vel,both device=/xs nxy=1,1 hann=15
echo "Convert sky-frequency to velocity "
echo "resolution of 1 km/s"
$<
# The line=velocity argument selects the range in terms of
# velocities in the data; in this case we take essentially
# the full range. The values used are determined by
# looking at the smauvspec plot from above, so you may run
# through this section a few times.
\rm -r $source.hco+vel.uv
uvaver vis=$source.hco+.uv line=velocity,390,-292,1 \
out=$source.hco+vel.uv
# Plot of Amplitude vs. channels and phase vs. channels
smauvspec vis=$source.hco+vel.uv \
interval=10000 options=nobase,avall \
axis=chan,both device=/xs nxy=1,1 hann=5
echo "Take out the continuum using linear fitting to"
echo "the line free channels"
\rm -r $source.hco+vel.uv.lin
# uvlin performs the continuum subtraction by fitting a polynomial
# to the real and imaginary parts of the line-free channels (here
# specified)
uvlin vis=$source.hco+vel.uv chans=1,100,290,390 \
out=$source.hco+vel.uv.lin\
order=1 mode=line
# chans=1,100,290,390 selects channels 1-100 and 290-390 for the
# continuum channels we use for subtraction -- from those we use
# linear interpolation to do the subtraction
# SgrA*: chans=5,80,150,220 (because there are different
# molecular features towards this source at various velocities than
# our target)
# 1E1740.7-2942: chans=1,100,290,390
echo "Press to continue."
$<
echo "Plot the continuum-free spectrum"
smauvspec vis=$source.hco+vel.uv.lin options=nobase,avall \
interval=1000 \
axis=vel,real device=/xs nxy=1,1 hann=15
echo "Press to continue."
$<
### CONSTRUCT LINE IMAGE CUBE ################################
# Invert the visibilities in the line of the source and
# deconvolve. Apply the model and look at the images.
##############################################################
echo "9. CONSTRUCT LINE IMAGE CUBE"
\rm -r $source.hco+.map $source.hco+.beam $source.hco+.icmp
# A note regarding inversion: see the manual page. The
# sidelobe suppression area (sup) is given in arcseconds, and
# gives the area around a source where inversion tries to
# suppress the sidelobes. The default shape is square, the
# default behavior is suppression across the field.
# The visibility weighting robustness parameter (robust=x)
# assigns weight to visibilities based on uv coverage.
# Natural weighting reduces the noise but degrades the
# beamshape.
# We can also change the fwhm parameter (fwhm=theta_maj,theta_min, pa);
# the SNR in the output is optimized when the parameter is set to the
# FWHM of features in the image. Careful: In deep cleaning, you
# can deconvolve noise this way too and make what looks like a convincing
# detection. One may set cutoff to stop the cleaning.
# The cutoff can be chosen based on the theoretical noise
# calculated from invert. Try cutoff= 2-3 sigma.
invert vis=$source.hco+vel.uv.lin map=$source.hco+.map \
beam=$source.hco+.beam \
imsize=512,512 cell=0.6 sup=0 \
line=vel,10,-150,2
# For S1, line=vel,10,-160,5, for S2 line=vel,10,-120,5 (we know
# where the features are from a previous CARMA observation)
\rm -r $source.hco+.icmp
clean map=$source.hco+.map beam=$source.hco+.beam \
out=$source.hco+.icmp gain=0.08 \
cutoff=.4 niters=5000 \
region=arcsec,'boxes(-25,-25,25,25)'
\rm -r $source.hco+.icln
restor model=$source.hco+.icmp beam=$source.hco+.beam \
map=$source.hco+.map \
out=$source.hco+.icln
cgdisp in=$source.hco+.icln type=c \
region=arcsec,'boxes(-25,-25,25,25)' \
xybin=1,1 device=/xs nxy=5,5 options=full,beambr,trlab,3val \
labtyp=arcsec,arcsec cols1=7 csize=0.5,0.5,0.5 \
slev=a,0.2 levs1=-5,3,4,8,16,32,64,128,256
echo "statistics on the channel images"
imstat in=$source.hco+.icln region=quarter \
axes=ra,dec
imstat in=$source.hco+.icln region='boxes(30,30,60,60)' \
axes=ra,dec
echo "Press to continue."
$<
cgdisp in=$source.hco+.icln type=pixel \
region=arcsec,'boxes(-25,-25,25,25)' \
xybin=1,1 device=/xs nxy=3,3 options=beambr,wedge,trlab,3val \
labtyp=arcsec,arcsec range=0,0,lin,2 cols1=7 csize=0.5,0.5,0.5
echo "Press to continue."
$<
# \rm -r $source.hco+.mom
# moment in=$source.hco+.icln region='image(3,6)' out=$source.hco+.mom axis=3 \
# mom=0
# cgdisp in=$source.hco+.mom device=/xs options=3val labtyp=arcsec
\rm -r $source.hco+.m0
moment in=$source.hco+.icln region='quarter(1,20)' \
out=$source.hco+.m0 axis=3 \
mom=0 clip=-1.2,1.2
smacgdisp in=$source.hco+.m0 type=pixel \
region=arcsec,'boxes(-30,-30,30,30)' \
xybin=1,1 device=/xs \
nxy=1,1 options=beambr,wedge,trlab,3val \
labtyp=arcsec,arcsec range=0,0,lin,2 cols1=7 csize=1,1,1
$<
\rm -r $source.hco+.m1
moment in=$source.hco+.icln region='quarter(1,20)' \
out=$source.hco+.m1 axis=3 \
mom=1 clip=-1.2,1.2
smacgdisp in=$source.hco+.m1 type=pixel \
region=arcsec,'boxes(-30,-30,30,30)' \
xybin=1,1 device=/xs \
nxy=1,1 options=beambr,wedge,trlab,3val \
labtyp=arcsec,arcsec range=-200,-100,lin,2 cols1=7 csize=1,1,1
# check the spectum from the image cube
# imspect in=$source.hco+.icln region=arcsec,'box(-3,-3,3,3)' device=/xs
echo " "
echo "The end"
echo " "
exit