Transient reduction pipeline

Transient Reduction Pipeline

2010 May 11 - T. Joseph W. Lazio

There are a class of observations for which only relatively simple data reduction steps are needed. One such example is that of transient observations, which are typically conducted in continuum mode and for which one is merely trying to determine a flux density of an essentially unresolved or only partially resolved source. This guide describes the steps in such a pipeline.

Warning: This guide was written at the time when the EVLA was still in its commissioning phase. As the instrument matures, specific steps taken here may need to be adjusted. Caveat emptor.

In order to process these data in a semi-automatic fashion, certain assumptions are made

• Data stored as an ALMA Science Data Model (ASDM) file on disk.
• Sources, listed in order of appearance, are structured as Option A or Option B
  • Option A: phase calibrator, target source, amplitude calibrator
  • Option B: amplitude calibrator, phase calibrator, target source

This guide should be set up so that a pipeline script can be extracted from it.

User Input

0. Initial set up, advice and warnings, and sanity checks.

# In CASA
import os.path, math
#
Version='v. 0.0 TJWL  2010-05-11'
CalModels='/home/casa/data/nrao/VLA/CalModels'
c=2.99792458E8
D=24.5
degrad=180/pi
arcsecond_deg=3600
#
print " "
print "EVLA Transient Reduction Pipeline"
print " "
print Version
print " "
print "Warning: This script was written at the time that the"
print "EVLA was undergoing commissioning.  Details may change."
print "Caveat emptor!"
print " "
print "Assumption: Data stored as an ALMA Science Data Model (ASDM)"
print "            file on disk."
print " "
print "Assumption: Sources, listed in order of appearance, are structured"
print "            as Option A or Option B"
print "            Option A: phase calibrator"
print "                      target source"
print "                      amplitude calibrator"
print " "
print "            Option B: amplitude calibrator"
print "                      phase calibrator"
print "                      target source"
print " "
#
print "Assumption: Calibrator model files stored in "
print CalModels
if not os.path.isdir(CalModels):
    raise IOError('Calibrator model files missing?  Stopping.')
#
print " "
acal_std_name = {}
acal_std_name['3C286'] = frozenset(['3C286','3C 286',
                                    '1328+30','1328+307','B1328+307','B1328+30',
                                    '1331+305','J1331+305','J1331+3030'])
acal_std_name['3C138'] = frozenset(['3C138','3C 138',
                                    '0518+16','0518+165','0518+1635','B0518+16','B0518+165','B0518+1635',
                                    '0521+166','0521+1638','J0521+166','J0521+1638'])
acal_std_name['3C147'] = frozenset(['3C147','3C 147',
                                    '0538+49','0538+498','0538+4949','B0538+49','B0538+498','B0538+4949',
                                    '0542+498','0542+4951','J0542+498','J0542+4951'])
acal_std_name['3C48'] = frozenset(['3C48','3C 48',
                                   '0134+32','0134+329','0134+3254','B0134+32','B0134+329','B0134+3254',
                                   '0137+331','0137+3309','J0137+331','J0137+3309'])

1. Obtain the name of the ASDM file, use it to construct the Measurement Set name.

# In CASA
ASDM_name=raw_input('ASDM file name:')
msname=ASDM_name+'.ms'

2. Determine the order in which sources appear in the observation.

# In CASA
obs_structure=raw_input('Structure of observations [A|B]:')
if obs_structure == 'A':
    sources = {'acal':'2', 'pcal':'0', 'target':'1'}
elif obs_structure == 'B':
    sources = {'acal':'0', 'pcal':'1', 'target':'2'}

3. Reference antenna. (Should be able to determine this from the measurement set!!)

# In CASA
refant=raw_input('Enter desired reference antenna [e.g., ea21]: ')

Processing Steps

1. Read the data from the ASDM file converting to a Measurement Set with importevla. Apply basic flagging operations (zeros, shadowing) here.

importevla(asdm=ASDM_name, vis=msname,
           flagzero=True, flagpol=True, shadow=True)

2. Flag first (dummy) scan.

flagdata(vis=msname,
         mode='manualflag',selectdata=T, scan='1')

3. Flag first 10 seconds of each scan.

flagdata(vis=msname,
         mode='quack',quackinterval=10,quackmode='beg',
         selectdata=F)

4. Extract various useful quantities from the data.

hdvalue=vishead(vis=msname, mode='get',
                hdkey='field', hdindex=sources['acal'])
acal_name=hdvalue[0]

spw_table=msname + '/SPECTRAL_WINDOW'
tb.open(spw_table)
freq_list=tb.getcol("REF_FREQUENCY")
channel_width=tb.getcol("CHAN_WIDTH")
num_chan=tb.getcol("NUM_CHAN")
tb.close(spw_table)
freq=freq_list[0]
if (1E9 < freq) and (freq < 2E9):
    band = 'L'
elif (2E9 < freq) and (freq < 4E9):
    band = 'S'
elif (4E9 < freq) and (freq < 8E9):
    band='C'
elif (8E9 < freq) and (freq < 12E9):
    band='X'
elif (12E9 < freq) and (freq < 40E9):
    band = 'K'
elif freq > 40E9:
    band = 'Q'
print "Observations are determined to be in the ", band, " band.\n"


wavelength=c/freq
primary_beam=(1.02*wavelength/D)*degrad   # HPBW, from Napier (1999)

tb.open(msname)
uvw=tb.getcol('UVW')
interval=tb.getcol('INTERVAL')
tb.close(msname)

b=[]
for u, v, w in zip(uvw[0], uvw[1], uvw[2]):
    b.append(math.sqrt(u*u+v*v)/wavelength)
bmax=max(b)
                                          # HPBW, synthesized beam,
                                          # from Bridle & Schwab (1999)
synthesized_beam=(1.2/bmax)*degrad*arcsecond_deg


##############################
# How much time-average
# smearing can be tolerated?
# Assume no more than time_smearing_loss peak
# intensity loss over half of
# primary beam.
# Follow Section 2 of Bridle
# & Schwab (1999).
#
time_smearing_loss=0.01
FoV=(primary_beam/math.sqrt(2))*arcsecond_deg
tau=math.sqrt(time_smearing_loss*1E9)*(synthesized_beam/FoV)

arbitrary_maximum=30
if tau > arbitrary_maximum:
    tau=arbitrary_maximum

dt=interval[0]
for dtprime in interval:
    if dtprime < dt:
        dt=dtprime

# tau is allowed value, dt is actual (minimum)
# (could be an issue if baseline-dependent
#  correlator accumulation used)
# make sure that tau is an integer
# multiple of dt
# 
tau=dt*math.floor(tau/dt)
if tau < dt:
    tau=dt

print "Data will be averaged in time."
print "Original time sampling [s]: ", dt
print "Averaging time [s]:         ", tau
print " "

##############################
# How much bandwidth
# smearing can be tolerated?
# Assume no more than band_smearing_loss peak
# intensity loss over half of
# primary beam.
# Follow Section 1 of Bridle
# & Schwab (1999).
# Assume square bandpass, no
# taper, expand resulting
# sine integral to lowest order
#
band_smearing_loss=0.01
eta_band=3.79
delta_nu=freq*(2/eta_band)*(synthesized_beam/FoV)*math.sqrt(18*band_smearing_loss)

# delta_nu is allowed value, channel_width
# is list of actual values
# make sure that delta_nu is an integer
# multiple of channel_width
# cheat

nchav=math.floor(delta_nu/min(channel_width[0]))
min_num_chan=max(1, num_chan)
if nchav > min_num_chan:
    nchav=min_num_chan

print "Data will be averaged in frequency."
print "Original channel width [kHz]: ", min(channel_width[0])/1E3
print "Averaged channel width [kHz]: ", (nchav*min(channel_width[0]))/1E3
print "Number of channels: ", nchav

5. Compress data in time.

   (You've checked that time average smearing isn't an issue ...?)

spw_table=msname + '/SPECTRAL_WINDOW'
tb.open(spw_table)
freq=tb.getcol("REF_FREQUENCY")[0]
tb.close(spw_table)
if (1E9 < freq) and (freq < 2E9):
    band = 'L'
elif (2E9 < freq) and (freq < 4E9):
    band = 'S'
elif (4E9 < freq) and (freq < 8E9):
    band='C'
elif (8E9 < freq) and (freq < 12E9):
    band='X'
elif (12E9 < freq) and (freq < 40E9):
    band = 'K'
elif freq > 40E9:
    band = 'Q'

outmsname=ASDM_name+'_split.ms'
split(vis=msname,outputvis=outmsname,
      datacolumn='data',
      timebin='10s')

6. Set flux density of amplitude calibrator.

  Need some heuristics here.
  Also could be an issue if the frequency setting is such that one should use a model image at a different band than one is observing, e.g., observing near the top of the C band where an X band model might be more appropriate.

for tstname in acal_std_name.keys():
    if acal_name in acal_std_name[tstname]:
        modimage=tstname
modimage=modimage + '_' + band + '.im'
setjy(vis=outmsname,field=sources['acal'],modimage=modimage)

7. Quick-n-dirty bandpass

Bcaltable=ASDM_name + '.B1'
bandpass(vis=outmsname,caltable=Bcaltable,
         field=sources['acal'],
         solint='inf',combine='scan',refant=refant,
         solnorm=T,
         bandtype='B',
         append=F)

8. amplitude and phase calibration

Gcaltable=ASDM_name + '.G1'
gaincal(vis=outmsname,caltable=Gcaltable,
        field=sources['acal'],
        solint='inf',refant=refant,
        minsnr=5,
        solnorm=F,
        gaintype='G',calmode='ap',
        append=F,
        gaintable=Bcaltable)

gaincal(vis=outmsname,caltable=Gcaltable,
        field=sources['pcal'],
        solint='inf',refant=refant,
        minsnr=5,
        solnorm=F,
        gaintype='G',calmode='ap',
        append=T,
        gaintable=Bcaltable)

9. Apply calibration.

fluxtable=ASDM_name + '.flux1'
fluxscale(vis=outmsname,caltable=Gcaltable,fluxtable=fluxtable,
          reference=sources['acal'],
          transfer=sources['pcal'],
          append=F)

applycal(vis=outmsname,
         field=sources['pcal'],
         gaintable=[Bcaltable, Gcaltable])

applycal(vis=outmsname,
         field=sources['target'],
         gaintable=[Bcaltable, Gcaltable])

10. Form target source measurement set.

vishead(vis=msname, mode='get',
        hdkey='field', hdindex=sources['target'])
target_name=vishead()

targetms=target_name + '.ms'
split(vis=outmsname,outputvis=targetms,
      datacolumn='corrected',
      field=sources['target'])

11. Image and CLEAN.

imagename=target_name + '.im'