# Calibration script for NGC 3256 SciVer ALMA data in CASA 4.7. # # Adopted/updated for CASA 4.7 from CASA ALMA Guide at # https://casaguides.nrao.edu/index.php/NGC3256_Band3_Calibration_for_CASA_4.3 # # Only final 'executive' CASA tasks are retained in the production-run script # - all the on-screen diagnostics and 'trial-error' steps are commented out. # # 2017/05/23, barta@asu.cas.cz, EU ARC, Czech node #--- Create a list of UIDs basename=['uid___A002_X1d54a1_X5','uid___A002_X1d54a1_X174', 'uid___A002_X1d54a1_X2e3','uid___A002_X1d5a20_X5', 'uid___A002_X1d5a20_X174','uid___A002_X1d5a20_X330'] #--- Print observation summary for name in basename: listobs(vis=name+'.ms') #--- Show antenna possitions plotants(vis=basename[0]+'.ms', figfile=basename[0]+'_plotants.png') #--- Create WVR and Tsys caltables. Remove the old ones before, for sure. for name in basename: os.system('rm -rf cal-tsys_'+name+'.calnew cal-*'+name+'.Wnew') # TSys gencal(vis=name+'.ms', caltype='tsys', caltable='cal-tsys_'+name+'.calnew') # WVR wvrgcal(vis=name+'.ms', caltable='cal-'+name+'.Wnew', toffset=-1, segsource=True, tie=["Titan,1037-295,NGC3256"], statsource="1037-295", smooth='2.88s') #--- Apriori flagging for name in basename: flagdata(vis=name+'.ms', flagbackup=False, mode = 'shadow') flagdata(vis=name+'.ms',mode='manual', autocorr=True) flagdata(vis=name+'.ms', mode='manual', flagbackup=False, intent='*POINTING*') flagdata(vis=name+'.ms', mode='manual', flagbackup=False, intent='*ATMOSPHERE*') # Keep the flagging done so far flagmanager(vis = name+'.ms', mode = 'save', versionname = 'Apriori') #--- Inspect the data - see large phase delays (wraps) for antenna DV07. # These phase distlrsions are Single Band Delays (SBD) for antenna DV07. # This holds just for first 3 UIDs taken on April 16 2011, the next day # the antenna was OK. #=== Commented out in production-run script ''' plotms(vis=basename[0]+'.ms', xaxis='channel', yaxis='phase', spw='3', antenna='PM03&DV07', correlation='XX', avgtime='1e8', showgui=True) ''' #--- Make correction calibration table for Apr 16 data for i in range(3): # loop over the first three ms's name=basename[i] os.system('rm -rf cal-'+name+'_del.K') # Caltable type SBD (single-band delays) gencal(vis=name+'.ms', caltable='cal-'+name+'_del.Knew', caltype='sbd', antenna='DV07', pol='X,Y', spw='1,3,5,7', parameter=[1.00, 1.10, -3.0, -3.0, -3.05, -3.05, -3.05, -3.05]) #--- Apply the corrections: SBD delays & WVR phase distorsion for Apr 16 data, # just the WVR corrections for Apr 17 data. for i in range(3): # loop over the first three data sets name=basename[i] applycal(vis=name+'.ms', flagbackup=False, spw='1,3,5,7', interp=['nearest','nearest'], gaintable=['cal-'+name+'_del.Knew', 'cal-'+name+'.Wnew']) for i in range(3,6): # loop over the last three data sets name=basename[i] applycal(vis=name+'.ms', flagbackup=False, spw='1,3,5,7', interp='nearest', gaintable='cal-'+name+'.Wnew') #--- See the corrections after calibration application: Switch between "DATA" # and "CORRECTED" column #=== Commented out in production-run script ''' plotms(vis=basename[0]+'.ms', xaxis='channel', yaxis='phase', spw='1', antenna='', correlation='XX', avgtime='1e8', coloraxis='baseline', avgscan=True, selectdata=True, field='1037*', showgui=True) ''' #--- Extract the corrected data to a new MS for name in basename: os.system('rm -rf '+name+'_K_WVR.ms*') split(vis=name+'.ms', outputvis=name+'_K_WVR.ms', datacolumn='corrected', spw='0~7') #--- Inspection of Tsys data (let it run as the output is to PNG files) for spw in ['1','3','5','7']: for name in basename: plotbandpass(caltable='cal-tsys_'+name+'.calnew', xaxis='freq', yaxis='amp', spw=spw, overlay='time', # also try overlay='antenna' plotrange=[0, 0, 40, 180], figfile='cal-tsys_per_spw_'+spw+'_'+name+'.png', interactive=False) #--- Scans 4,5 and 9 at antenna DV04 in SPW7 are affected, let us flag them flagdata(vis='uid___A002_X1d54a1_X174_K_WVR.ms', mode='manual', antenna='DV04', flagbackup=True, scan='4,5,9', spw='7') #--- Apply the TSys calibration, extra for each field for name in basename: for field in ['Titan','1037*','NGC*']: applycal(vis=name+'_K_WVR.ms', spw='1,3,5,7', flagbackup=False, field=field, gainfield=field, interp='linear,spline', gaintable=['cal-tsys_'+name+'.calnew']) #--- Extract corrected science data to separate MS - leave the not usefull # WVR and TSys calibration data behind for name in basename: os.system('rm -rf '+name+'_line.ms*') split(vis=name+'_K_WVR.ms', outputvis=name+'_line.ms', datacolumn='corrected', spw='1,3,5,7') #--- Concatenate the data in the one big MS # Prepare the list of corrected partial MSs comvis=[] for name in basename: comvis.append(name+'_line.ms') # Just for sure remove (if it exists) os.system('rm -rf ngc3256_line.ms*') # Make a new cocatenated MS concat(vis=comvis, concatvis='ngc3256_line.ms') #--- Data inspection listobs('ngc3256_line.ms') # Show spectrum - see the edge channels #=== Commented out in production-run script ''' plotms(vis='ngc3256_line.ms', xaxis='channel', yaxis='amp', averagedata=True, avgbaseline=True, avgtime='1e8', avgscan=True, showgui=True) ''' #--- Remove the edge channels flagdata(vis='ngc3256_line.ms', flagbackup=True, spw='*:0~16,*:125~127') #=== Commented out in production-run script # Titan has an increase of intesity becoause blending with Saturn rings # on the next day... ''' plotms(vis='ngc3256_line.ms', xaxis='time', yaxis='amp', averagedata=True, avgchannel='128', coloraxis='field', iteraxis='spw', showgui=True) ''' # ...this can be seen here by source size increase. ''' plotms(vis='ngc3256_line.ms', xaxis='uvdist', yaxis='amp', field='Titan', averagedata=T, avgchannel='128', avgtime='1e8', coloraxis='scan') ''' #=== /End of diagnostic plots #--- Flag the Saturn-ring blended data flagdata(vis = 'ngc3256_line.ms', flagbackup=True, timerange='>2011/04/16/12:00:00', field='Titan') #--- Give the Titan an ephemeris - not needed anymore as new ASDM format keeps # it internally fixplanets(vis='ngc3256_line.ms', field='Titan', fixuvw=True) #--- Further data inspection plotms(): # * Baselines with DV07 have very high amplitudes in spw 3, correlation YY # * Baselines with DV08 have very low amplitudes in spw 3, correlation YY, but # only for the last observation # * Baselines with PM03 have low amplitudes at 2011/04/17/02:15:00 for spw 0 # * Baselines with PM03 have low amplitudes at 2011/04/16/04:15:15 for spw 2 # and 3 # * For the second day, the baseline DV10-PM03 is corrupted # # We flag the corresponding data flagdata(vis='ngc3256_line.ms', flagbackup=True, spw='3', correlation='YY', mode='manual', antenna='DV07', timerange='') flagdata(vis='ngc3256_line.ms', flagbackup=True, spw='3', correlation='YY', mode='manual', antenna='DV08', timerange='>2011/04/17/03:00:00') flagdata(vis='ngc3256_line.ms', flagbackup=True, spw='0', correlation='', mode='manual', antenna='PM03', timerange='2011/04/17/02:15:00~02:15:50') flagdata(vis='ngc3256_line.ms', flagbackup=True, spw='2,3', correlation='', mode='manual', antenna='PM03', timerange='2011/04/16/04:13:50~04:18:00') flagdata(vis='ngc3256_line.ms', flagbackup=True, spw='', correlation='', mode='manual', antenna='PM03&DV10', timerange='>2011/04/16/15:00:00') #--- Inspect bandpass (B) calibrator #=== Commented out in production-run script ''' plotms(vis='ngc3256_line.ms', xaxis='freq', yaxis='phase', selectdata=True, field='1037*', avgtime='1E6', avgscan=True, coloraxis='baseline', iteraxis='antenna', showgui=True) plotms(vis='ngc3256_line.ms', xaxis='time', yaxis='phase', selectdata=True, field='1037*', spw='0:30~90', avgchannel='128', avgscan=True, coloraxis='baseline', iteraxis='antenna', showgui=True) ''' #+++++++++++++ Inspection of primary phase (P) calibrator (1037-295) after # application of calibration tables shows outliers fo some scans and particular # antennas. With this apriori knowledge (found by trial-error, of course) we # can flag the data right now to avoid repeated calibration steps. #--- Flag corrupted data flagdata(vis='ngc3256_line.ms', mode='manual', timerange='2011/04/16/04:13:35~04:13:45', flagbackup=True) # tget(): # copy-paste last args of flagdata() call, we will change just some of them tget(flagdata) timerange='2011/04/16/05:21:13~05:21:19' flagdata() timerange='2011/04/16/04:16:40~04:16:49' flagdata() timerange='2011/04/16/04:14:00~04:17:10'; antenna='PM03' flagdata() timerange='2011/04/17/01:52:20~01:53:10'; antenna='DV10' flagdata() timerange='2011/04/17/00:35:30~01:20:20'; antenna='DV04'; spw='3' flagdata() #--- Short time-scale balancing of phase variations for the B calibrator gaincal(vis='ngc3256_line.ms', caltable='cal-ngc3256.G1', spw='*:40~80', field='1037*', selectdata=True, solint='int', refant='DV04', calmode='p') #--- Inspect the phase caltables for B calibrator #=== Commented out in production-run script ''' plotcal(caltable = 'cal-ngc3256.G1', xaxis = 'time', yaxis = 'phase', poln='X', plotsymbol='o', plotrange = [0,0,-180,180], iteration = 'spw', figfile='cal-phase_vs_time_XX.G1.png', subplot = 221, showgui=True) plotcal(caltable = 'cal-ngc3256.G1', xaxis = 'time', yaxis = 'phase', poln='Y', plotsymbol='o', plotrange = [0,0,-180,180], iteration = 'spw', figfile='cal-phase_vs_time_YY.G1.png', subplot = 221, showgui=True) ''' #--- Do bandpass calibration ('spectral balancing') separatelly for Apr 16 # and 17 data. Make phase corrections for B calibrator 'on the fly' # (instead of calling applycal() beforehand) - use parameter 'gaintable' # Apr 16 data bandpass(vis = 'ngc3256_line.ms', caltable = 'cal-ngc3256.B1', gaintable = 'cal-ngc3256.G1', timerange='<2011/04/16/15:00:00', field = '1037*', minblperant=3, minsnr=2, solint='inf', combine='scan,obs', bandtype='B', fillgaps=1, refant = 'DV04', solnorm=True) # Apr 17 data - append to the caltable made in previous step bandpass(vis = 'ngc3256_line.ms', caltable = 'cal-ngc3256.B1', gaintable = 'cal-ngc3256.G1', timerange='>2011/04/16/15:00:00', field = '1037*', minblperant=3, minsnr=2, solint='inf', combine='scan,obs', bandtype='B', fillgaps=1, refant = 'DV04', solnorm=True, append=True) #--- Plot bandpass solutions: Phase & Amplitude #=== Commented out in production-run script ''' plotbandpass(caltable = 'cal-ngc3256.B1', xaxis='freq', yaxis='phase', plotrange = [0,0,-70,70], overlay='antenna', figfile='bandpass.B1.png') plotbandpass(caltable = 'cal-ngc3256.B1', xaxis='freq', yaxis='amp', overlay='antenna', figfile='bandpass.B2.png') ''' #--- Absolute flux calibration for Titan setjy(vis='ngc3256_line.ms', field='Titan', standard='Butler-JPL-Horizons 2012', spw='0,1,2,3', usescratch=False) #--- Long time-scale (complex) gain calibration for primary and secondary # phase (P) calibrators - 1037-295 and Titan: make table G2 # Use the B-calibration on the fly (use 'gaintable' parameter). gaincal(vis = 'ngc3256_line.ms', caltable = 'cal-ngc3256.G2', spw = '*:16~112', field = '1037*,Titan', minsnr=1.0, solint= 'inf', selectdata=True, solnorm=False, refant = 'DV04', gaintable = 'cal-ngc3256.B1', calmode = 'ap') #--- Inspect the G2 calibration table #=== Commented out in production-run script ''' plotcal(caltable = 'cal-ngc3256.G2', xaxis = 'time', yaxis = 'phase', poln='X', plotsymbol='o', plotrange = [0,0,-180,180], iteration= 'spw', figfile='cal-phase_vs_time_XX.G2.png', subplot = 221, showgui=True) plotcal(caltable = 'cal-ngc3256.G2', xaxis = 'time', yaxis = 'phase', poln='Y', plotsymbol='o', plotrange = [0,0,-180,180], iteration= 'spw', figfile='cal-phase_vs_time_YY.G2.png', subplot = 221, showgui=True) plotcal(caltable = 'cal-ngc3256.G2', xaxis = 'time', yaxis = 'amp', poln='X', plotsymbol='o', plotrange = [], iteration = 'spw', figfile='cal-amp_vs_time_XX.G2.png', subplot = 221, showgui=True) plotcal(caltable = 'cal-ngc3256.G2', xaxis = 'time', yaxis = 'amp', poln='Y', plotsymbol='o', plotrange = [], iteration = 'spw', figfile='cal-amp_vs_time_YY.G2.png', subplot = 221, showgui=True) ''' #--- Bootstrapping of the absolute flux density for the phase calibrator # 1037-295 from Titan fluxscale( vis="ngc3256_line.ms", caltable="cal-ngc3256.G2", fluxtable="cal-ngc3256.G2.flux", reference="Titan", transfer="1037*", refspwmap=[0,1,1,1]) #--- Apply the calibrations to the science target and the phase calibrator applycal( vis='ngc3256_line.ms', flagbackup=False, field='NGC*,1037*', interp=['nearest','nearest'], gainfield = ['1037*', '1037*'], gaintable=['cal-ngc3256.G2.flux', 'cal-ngc3256.B1']) #--- Inspect calibrated data: Show the spectrum of science target #=== Commented out in production-run script ''' plotms(vis='ngc3256_line.ms', xaxis='freq', yaxis='amp', ydatacolumn='corrected', selectdata=True, field='NGC*', averagedata=True, avgchannel='', avgtime='10000', avgscan=True, avgbaseline=True, coloraxis='spw', showgui=True) ''' #--- Inspect primary P calibrator #=== Commented out in production-run script # amplitude vs time ''' plotms(vis='ngc3256_line.ms', xaxis='time', yaxis='amp', ydatacolumn='corrected', selectdata=True, field='1037*', averagedata=True, avgchannel='128', avgtime='', avgscan=False, avgbaseline=False, coloraxis='spw', showgui=True) ''' # phase vs time ''' plotms(vis='ngc3256_line.ms', xaxis='freq', yaxis='phase', ydatacolumn='corrected', selectdata=True, field='1037*', averagedata=True, avgchannel='', avgtime='1e8', avgscan=True, avgbaseline=False, coloraxis='baseline', showgui=True) ''' #++++++++++ The above inspection of P calibrator would find the corrupted data # unless we flagged them in advance (what we did). # The following sequence - repetition of calibration steps with late flagged # data is thus unnecessary now and is commented out: ''' #--- Flag corrupted data flagdata(vis='ngc3256_line.ms', mode='manual', timerange='2011/04/16/04:13:35~04:13:45', flagbackup=True) # copy-paste last args of flagdata() call, we will change just some of them tget(flagdata) timerange='2011/04/16/05:21:13~05:21:19' flagdata() timerange='2011/04/16/04:16:40~04:16:49' flagdata() timerange='2011/04/16/04:14:00~04:17:10'; antenna='PM03' flagdata() timerange='2011/04/17/01:52:20~01:53:10'; antenna='DV10' flagdata() timerange='2011/04/17/00:35:30~01:20:20'; antenna='DV04'; spw='3' flagdata() #--- Re-initialise the "MODEL" column for calibrators after flagging delmod('ngc3256_line.ms') #--- Repeat all the calibration with data shrinked by flagging gaincal(vis='ngc3256_line.ms', caltable='cal-ngc3256.G1n', spw='*:40~80', field='1037*', selectdata=True, solint='int', refant='DV04', calmode='p') bandpass(vis = 'ngc3256_line.ms', caltable = 'cal-ngc3256.B1n', gaintable = 'cal-ngc3256.G1n', timerange='<2011/04/16/15:00:00', field = '1037*', minblperant=3, minsnr=2, solint='inf', combine='scan,obs', bandtype='B', fillgaps=1, refant = 'DV04', solnorm=True) bandpass(vis = 'ngc3256_line.ms', caltable = 'cal-ngc3256.B1n', gaintable = 'cal-ngc3256.G1n', timerange='>2011/04/16/15:00:00', field = '1037*', minblperant=3, minsnr=2, solint='inf', combine='scan,obs', bandtype='B', fillgaps=1, refant = 'DV04', solnorm=True, append=True) setjy(vis='ngc3256_line.ms', field='Titan', standard='Butler-JPL-Horizons 2012', spw='0,1,2,3',usescratch=False) gaincal(vis = 'ngc3256_line.ms', caltable = 'cal-ngc3256.G2n', spw = '*:16~112', field = '1037*,Titan', minsnr=1.0, solint= 'inf', selectdata=True, solnorm=False, refant = 'DV04', gaintable = 'cal-ngc3256.B1n', calmode = 'ap') fluxscale(vis="ngc3256_line.ms", caltable="cal-ngc3256.G2n", fluxtable="cal-ngc3256.G2n.flux", reference="Titan", transfer="1037*", refspwmap=[0,1,1,1]) applycal(vis='ngc3256_line.ms', flagbackup=False, field='NGC*,1037*', interp=['nearest','nearest'], gainfield = ['1037*', '1037*'], gaintable=['cal-ngc3256.G2n.flux', 'cal-ngc3256.B1n']) #--- Inspect the phase calibrator again plotms(vis='ngc3256_line.ms', xaxis='time', yaxis='amp', ydatacolumn='corrected', selectdata=True, field='1037*', averagedata=True, avgchannel='128', avgtime='', avgscan=False, avgbaseline=False, coloraxis='spw', showgui=True) plotms(vis='ngc3256_line.ms', xaxis='freq', yaxis='phase', ydatacolumn='corrected', selectdata=True, field='1037*', averagedata=True, avgchannel='', avgtime='1e8', avgscan=True, avgbaseline=False, coloraxis='baseline', showgui=True) # now OK #--- /end of repeted steps ''' #--- Imaging of calibrators os.system('rm -rf result-phasecal_cont*') # Primary P-calibrator - Inverse FFT and deconvolution in-one clean(vis='ngc3256_line.ms', imagename='result-phasecal_cont', field='1037*', spw='*:20~120', selectdata=True, mode='mfs', niter=500, gain=0.1, threshold='0.75mJy', psfmode='hogbom', interactive=False, mask=[62, 62, 67, 67], imsize=128, cell='1arcsec', weighting='briggs', robust=0.0, nterms=2) #--- make image statistics calstat=imstat(imagename='result-phasecal_cont.image.tt0', region='', box='85,8,120,120') rms=(calstat['rms'][0]) print '>> rms in phase calibrator image: '+str(rms) calstat=imstat(imagename='result-phasecal_cont.image.tt0', region='') peak=(calstat['max'][0]) print '>> Peak in phase calibrator image: '+str(peak) print '>> Dynamic range in phase calibrator image: '+str(peak/rms) #--- View the primary P-calibrator: set output to PNG file imview(raster={'file': 'result-phasecal_cont.image.tt0', 'colorwedge':T, 'range':[-0.004, 0.250], 'scaling':-1.5, 'colormap':'Rainbow 2'}, out='result-phasecal_map.png', zoom=1) #--- Apply calibrations for Titan applycal(vis='ngc3256_line.ms', flagbackup=False, field='Titan', interp=['nearest', 'nearest'], gainfield = ['Titan', '1037*'], gaintable=['cal-ngc3256.G2.flux', 'cal-ngc3256.B1']) os.system('rm -rf result-ampcal_cont*') #--- Titan - Inverse FFT and deconvolution in-one clean(vis='ngc3256_line.ms', imagename='result-ampcal_cont', field='Titan', spw='0:20~120,1:20~120', mode='mfs', niter=200, threshold='5mJy', psfmode='hogbom', mask=[62, 62, 67, 67], imsize=128, cell='1arcsec', weighting='briggs', robust=0.0, interactive=False) #--- make image statistics calstat=imstat(imagename="result-ampcal_cont.image", region="",box="85,8,120,120") rms=(calstat['rms'][0]) print ">> rms in amp calibrator image: "+str(rms) calstat=imstat(imagename="result-ampcal_cont.image",region="") peak=(calstat['max'][0]) print ">> Peak in amp calibrator image: "+str(peak) print ">> Dynamic range in amp calibrator image: "+str(peak/rms) #--- View Titan: set output to PNG file imview(raster={'file': 'result-ampcal_cont.image', 'colorwedge':T, 'range':[-0.02, 0.250], 'scaling':-1.5, 'colormap':'Rainbow 2'}, out='result-ampcal_map.png', zoom=1) os.system('rm -rf ngc3256_line_target.ms*') #--- Extract just the science target NGC 3256 into separate MS split(vis='ngc3256_line.ms', outputvis='ngc3256_line_target.ms', field='NGC*')