# ALMA Data Reduction Script

# Self-calibration for uid___A002_Xc89480_X10c4c.ms target RXCJ2341.1+0018
# RXCJ2341.split contains the target only, after all pipeline calibration.
# The line is weak so the line spw has been averaged to the same resoluton as the continuum
# and a separate file contains the central quarter of the line spw with less averaging

contvis='RXCJ0439.split'
linevis='RXCJ0439_line.split'
targ ='RXCJ0439.0+0520'

antref='DA47'

# imaging parameters
active=True            # whether to clean interactively or not)

c = '0.05arcsec'       # pixel size  
isz = 300              # image size in pixels
spmode='mfs'           # continuum mode
contchans='0,1,2,3:0~44;70~119'     #line-free channels
r = 0.5                # robust parameter value  
sig = '20.e-6Jy'      # expected stopping threshold
mymask='user'

# calibration parameters
p1solint='60s'
p1mode='p'
p1gtype='G'
p1apply='calflag'

a1solint='inf'
a1mode='a'
a1gtype='T'
a1apply='calonly'
ap1gtables=[contvis+'.p1',contvis+'.a1'] # any gaintables applied in gaincal plus output caltable

# line data parameters
linecontchans='0:0~399;560~959'
linespmode='cube'
a=400           # start
b=2             # averaging chans 
n=80            # nchans
mylinemask='auto-multithresh'
linesig='0.2mJy'
liner=1.5
nthresh = 4.    # rms threshold for automasking
sthresh=2.      # sidelobe rejection threshold 

# box sizes for measuring rms and peak
rbox='10,10'+','+str(isz-10)+','+str(isz*0.2)
pbox= str(isz*0.2)+','+str(isz*0.2)+','+str(isz*0.8)+','+str(isz*0.8)

thesteps = []
step_title = {0: 'listobs, plotms',
              1: 'flagmanager, optional flagging',
              2: 'First continuum image',
              3: 'First round of phase self-cal',
              4: 'Apply solutions and re-image',
              5: 'First round of amp self-cal',
              6: 'Apply solutions and re-image',
              7: 'Apply solutions to line and plotms',
              8: 'Subtract continuum',
              9: 'Make line cube'
}

try:
  print 'List of steps to be executed ...', mysteps
  thesteps = mysteps
except:
  print 'global variable mysteps not set.'
if (thesteps==[]):
  thesteps = range(0,len(step_title))
  print 'Executing all steps: ', thesteps

# The Python variable 'mysteps' will control which steps
# are executed when you start the script using
#   execfile('scriptForCalibration.py')
# e.g. setting
#   mysteps = [2,3,4]
# before starting the script will make the script execute
# only steps 2, 3, and 4
# Setting mysteps = [] will make it execute all steps.

mystep = 0
if(mystep in thesteps):
  casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
  print 'Step ', mystep, step_title[mystep]

  clearcal(vis=contvis)   # remove any previous self-calibration
  clearcal(vis=linevis)

  listobs(vis=contvis,
          overwrite=True,
          listfile=contvis+'.listobs')

  listobs(vis=linevis,
          overwrite=True,
          listfile=linevis+'.listobs')

# This plotms is to identify line channels to exclude from continuum
  for s in ['0','1','2','3']:
    plotms(vis=contvis,
           xaxis='channel',
           yaxis='amp',
           spw=s,
           avgtime='3600', 
           avgscan=True,   # average all times
           avgbaseline=True,
           overwrite=True,
           plotfile=contvis+'spw'+s+'amp-chan.png')

  plotms(vis=contvis,
           xaxis='uvdist',
           yaxis='amp',
           avgchannel='128',
           avgtime='3600', # average per scan
           overwrite=True,
           plotfile=contvis+'_amp-uvdist.png')

# This plotm is for you to inspect the phase v. baseline to help choose solints
# investigate and zoom interactively
  plotms(vis=contvis,
           xaxis='time',
           yaxis='phase',
           spw='1',
           antenna=antref+'&*',
           avgchannel='128',
           coloraxis='corr',
           iteraxis='baseline')  

mystep = 1
if(mystep in thesteps):
  casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
  print 'Step ', mystep, step_title[mystep]

  flagmanager(vis=contvis,
              mode='save',
              versionname='post-split')

  flagmanager(vis=linevis,
              mode='save',
              versionname='post-split')
# add flagging commands if wanted

mystep = 2
if(mystep in thesteps):
  casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
  print 'Step ', mystep, step_title[mystep]

  os.system('rm -rf '+targ+'.clean*')
  tclean(vis = contvis,
         imagename = targ+'.clean',
         specmode = spmode,
         spw=contchans,
         imsize = [isz, isz],
         cell = c,
         niter = 1000,           
         weighting = 'briggs',
         robust = r,
         mask = '',
         usemask=mymask,
         gridder = 'standard',
         pbcor = True,
         threshold = sig,
         interactive = active
         )

  rms=imstat(imagename=targ+'.clean.image', box=rbox)['rms'][0]
  peak=imstat(imagename=targ+'.clean.image', box=pbox)['max'][0]
  print targ+'.clean.image: Peak %7.3f mJy/bm, rms  %6.3f mJy, S/N %5.1f' % (peak*1000., rms*1000., peak/rms)

  ft(vis = contvis,
     model=targ+'.clean.model',
     nterms=1,
     usescratch=True)

mystep = 3
if(mystep in thesteps):
  casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
  print 'Step ', mystep, step_title[mystep]

  os.system('rm -rf '+contvis+'.p1')
  gaincal(vis = contvis,
          caltable=contvis+'.p1',
          solint=p1solint,
          refant=antref,
          calmode=p1mode,
          minsnr=2,           # relax to avoid failing solutions if previous calibration or model, not data are bad
          gaintype=p1gtype)

  plotcal(caltable=contvis+'.p1',
          xaxis='time',
          yaxis='phase',
          subplot=551,
          iteration='antenna',
          antenna='DA*',
          figfile=contvis+'.p1_DA.png')

  plotcal(caltable=contvis+'.p1',
          xaxis='time',
          yaxis='phase',
          subplot=551,
          iteration='antenna',
          antenna='DV*,PM*',
          figfile=contvis+'.p1_DVPM.png')

mystep = 4
if(mystep in thesteps):
  casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
  print 'Step ', mystep, step_title[mystep]

  flagmanager(vis = contvis,
              mode='save',
              versionname='beforeapplycal')

  applycal(vis = contvis,
           gaintable=contvis+'.p1',
           calwt=False,
           applymode=p1apply,
           flagbackup=False)

  os.system('rm -rf '+targ+'.cleanp1*')
  tclean(vis = contvis,
         imagename = targ+'.cleanp1',
         specmode = spmode,
         deconvolver='mtmfs',     # if we are going to do amp self-cal
         nterms=2,                # fit a linear slope
         spw=contchans,
         imsize = [isz, isz],
         cell = c,
         niter = 1000,           
         weighting = 'briggs',
         robust = r,
         mask = '',
         usemask=mymask,
         gridder = 'standard',
         pbcor = True,
         threshold = sig,
         interactive = active
         )

  rms=imstat(imagename=targ+'.cleanp1.image.tt0', box=rbox)['rms'][0]
  peak=imstat(imagename=targ+'.cleanp1.image.tt0', box=pbox)['max'][0]
  print targ+'.cleanp1.image: Peak %7.3f mJy/bm, rms  %6.3f mJy, S/N %5.1f' % (peak*1000., rms*1000., peak/rms)

  ft(vis = contvis,
     model=[targ+'.cleanp1.model.tt0',targ+'.cleanp1.model.tt1'],
     nterms=2,
     usescratch=True)

mystep = 5
if(mystep in thesteps):
  casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
  print 'Step ', mystep, step_title[mystep]

  os.system('rm -rf '+contvis+'.a1')
  gaincal(vis = contvis,
          caltable=contvis+'.a1',
          solint=a1solint,
          refant=antref,
          calmode=a1mode,
          gaintable=contvis+'.p1',   # existing phase solutions
          gaintype=a1gtype)

  plotcal(caltable=contvis+'.a1',
          xaxis='time',
          yaxis='amp',
          subplot=551,
          iteration='antenna',
          antenna='DA*',
          figfile=contvis+'.a1_DA.png')

  plotcal(caltable=contvis+'.a1',
          xaxis='time',
          yaxis='amp',
          subplot=551,
          iteration='antenna',
          antenna='DV*,PM*',
          figfile=contvis+'.a1_DVPM.png')

mystep = 6
if(mystep in thesteps):
  casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
  print 'Step ', mystep, step_title[mystep]

  applycal(vis = contvis,
           gaintable=ap1gtables,   
           calwt=False,
           applymode=p1apply,
           flagbackup=False)

  os.system('rm -rf '+targ+'.cleanap1*')
  tclean(vis = contvis,
         imagename = targ+'.cleanap1',
         specmode = spmode,
         deconvolver='mtmfs',     # if we are going to do amp self-cal
         nterms=2,                # fit a linear slope
         spw=contchans,
         imsize = [isz, isz],
         cell = c,
         niter = 1000,           
         weighting = 'briggs',
         robust = r,
         usemask=mymask,
         mask = '',
         gridder = 'standard',
         pbcor = True,
         threshold = sig,
         interactive = active
         )

  rms=imstat(imagename=targ+'.cleanap1.image.tt0', box=rbox)['rms'][0]
  peak=imstat(imagename=targ+'.cleanap1.image.tt0', box=pbox)['max'][0]
  print targ+'.cleanap1.image: Peak %7.3f mJy/bm, rms  %6.3f mJy, S/N %5.1f' % (peak*1000., rms*1000., peak/rms)

  ft(vis = contvis,
     model=[targ+'.cleanap1.model.tt0',targ+'.cleanap1.model.tt1'],
     nterms=2,
     usescratch=True)

mystep = 7
if(mystep in thesteps):
  casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
  print 'Step ', mystep, step_title[mystep]

  applycal(vis = linevis,
           gaintable=ap1gtables,   
           calwt=False,
           applymode=p1apply,
           flagbackup=False)

  plotms(vis=linevis,
           ydatacolumn='corrected',
           xaxis='channel',
           yaxis='amp',
           uvrange='<200',
           avgtime='3600', 
           avgscan=True, 
           avgbaseline=True, 
           overwrite=True,
           plotfile=linevis+'_amp-channel.png')


mystep = 8
if(mystep in thesteps):
  casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
  print 'Step ', mystep, step_title[mystep]

  os.system('rm -rf '+linevis+'.contsub')
  mstransform(vis = linevis,
              outputvis=linevis+'.contsub',
              douvcontsub=True,
              fitorder=1,
              fitspw=linecontchans)

mystep = 9
if(mystep in thesteps):
  casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
  print 'Step ', mystep, step_title[mystep]

# don't forget to select all channels if you want when masking
  os.system('rm -rf '+targ+'_line.clean*')
  tclean(vis = linevis+'.contsub',
         imagename = targ+'_line.clean',
         specmode = linespmode,
         start=a,
         width=b,
         nchan=n,
         imsize = [isz, isz],
         cell = c,
         niter = 2000,           
         weighting = 'briggs',
         robust = liner,
         mask = '',
         usemask=mylinemask,
         noisethreshold=nthresh,
         negativethreshold=nthresh,
         sidelobethreshold=sthresh,
         gridder = 'standard',
         pbcor = True,
         threshold = linesig,
         interactive = active
         )