###################################################################### # # # Reduction script for 3C129, 6cm VLA # # Adapted by David Whysong from Steve Myers' Jupiter script # # Updated DHW 2007-11-14 # # Last Updated GvM 2007-12-04 minor corrections # # # ###################################################################### import time import os # #===================================================================== # # This script has some interactive commands: scriptmode = True # if you are running it and want it to stop during interactive parts. scriptmode = True #===================================================================== # # Set up some useful variables - these will be set during the script # also, but if you want to restart the script in the middle here # they are in one place: pathname=os.environ.get('AIPSPATH').split()[0] prefix='3c129_6cm' msfile = prefix + '.ms' gtable = prefix + '.gcal' ftable = prefix + '.fluxscale' atable = prefix + '.accum' srcsplitms = prefix + '.split.ms' imname1 = prefix + '.clean1' clnimage1 = imname1+'.image' clnmodel1 = imname1+'.model' clnresid1 = imname1+'.residual' clnmask1 = imname1+'.clean_interactive.mask' selfcaltab1 = srcsplitms + '.selfcal1' imname2 = prefix + '.clean2' clnimage2 = imname2+'.image' clnmodel2 = imname2+'.model' clnresid2 = imname2+'.residual' clnmask2 = imname2+'.clean_interactive.mask' selfcaltab2 = srcsplitms + '.selfcal2' smoothcaltab2 = srcsplitms + '.smoothcal2' imname3 = prefix + '.clean3' clnimage3 = imname3+'.image' clnmodel3 = imname3+'.model' clnresid3 = imname3+'.residual' clnmask3 = imname3+'.clean_interactive.mask' # #===================================================================== # # Get to path to the CASA home and stip off the name pathname=os.environ.get('AIPSPATH').split()[0] # This is where the UVFITS data will be #fitsdata=pathname+'/data/nrao/VLA/planets_6cm.fits' #fitsdata='/home/sandrock2/smyers/NAUG2/Data/VLA_CONT/FLUX99-6CM.CBAND' #fitsdata='3c129_6cm.uvfits' fitsdata='planets_6cm.fits' # The prefix to use for all output files prefix='3c129_6cm' # Clean up old files os.system('rm -rf '+prefix+'*') # #===================================================================== # Data Import and List #===================================================================== # # Import the data from FITS to MS # print '--Import--' # Safest to start from task defaults default('importvla') print "Use importvla to read VLA archive files and make an MS" # Set up the MS filename and save as new global variable msfile = prefix + '.ms' print "MS will be called "+msfile importvla(vis=msfile, bandname='C', archivefiles=['AT166_1', 'AT166_2', 'AT166_3']) #===================================================================== # # List a summary of the MS # print '--Listobs--' # Don't default this one and make use of the previous setting of # vis. Remember, the variables are GLOBAL! print "Use listobs to print verbose summary to logger" # You may wish to see more detailed information, in this case # use the verbose = True option listobs(verbose=T) # You should get in your logger window and in the casapy.log file # something like: # # Wed Nov 28 21:39:11 2007 INFO listobs::::: # ############################################### # # Wed Nov 28 21:39:11 2007 INFO listobs::::: # ### Begin Task: listobs ### # # Wed Nov 28 21:39:12 2007 INFO listobs::ms::summary: # MeasurementSet Name: /users/dwhysong/casa/casa-test/3c129_6cm.ms MS Version 2 # # Observer: unavailable Project: AT166 # Observation: VLA # # Telescope Observation Date Observer Project # VLA [ 4.28184e+09, 4.28184e+09]unavailable AT166 # VLA [ 4.28184e+09, 4.28187e+09]unavailable AT166 # VLA [ 4.29055e+09, 4.29057e+09]unavailable AT166 # # Wed Nov 28 21:39:12 2007 INFO listobs::ms::summary: # Data records: 677430 Total integration time = 8.7342e+06 seconds # Observed from 07:41:45 to 09:51:45 # # Wed Nov 28 21:39:12 2007 INFO listobs::ms::summary: # ObservationID = 0 ArrayID = 0 # Date Timerange Scan FldId FieldName SpwIds # 25-Jul-1994/07:41:45.0 - 07:43:25.0 1 0 2345-167 [0, 1] # # Wed Nov 28 21:39:13 2007 INFO listobs::ms::summary: # ObservationID = 1 ArrayID = 0 # Date Timerange Scan FldId FieldName SpwIds # 25-Jul-1994/08:00:25.0 - 08:02:35.0 2 1 A2597 [0, 1] # 08:21:05.0 - 08:22:45.0 3 0 2345-167 [0, 1] # 08:39:45.0 - 08:41:55.0 4 1 A2597 [0, 1] # 09:00:25.0 - 09:02:05.0 5 0 2345-167 [0, 1] # 09:19:05.0 - 09:21:15.0 6 1 A2597 [0, 1] # 09:39:45.0 - 09:41:25.0 7 0 2345-167 [0, 1] # 09:58:25.0 - 10:00:35.0 8 1 A2597 [0, 1] # 10:19:05.0 - 10:20:45.0 9 0 2345-167 [0, 1] # 10:37:45.0 - 10:39:55.0 10 1 A2597 [0, 1] # 10:58:25.0 - 11:00:05.0 11 0 2345-167 [0, 1] # 11:17:05.0 - 11:19:15.0 12 1 A2597 [0, 1] # 11:37:45.0 - 11:39:25.0 13 0 2345-167 [0, 1] # 11:48:25.0 - 11:50:05.0 14 2 0420+417 [0, 1] # 12:12:55.0 - 12:23:05.0 15 3 3C129 [0, 1] # 12:27:25.0 - 12:29:05.0 16 2 0420+417 [0, 1] # 12:50:05.0 - 13:00:15.0 17 3 3C129 [0, 1] # 13:04:35.0 - 13:06:15.0 18 2 0420+417 [0, 1] # 13:29:05.0 - 13:39:15.0 19 3 3C129 [0, 1] # 13:43:35.0 - 13:45:15.0 20 2 0420+417 [0, 1] # # Wed Nov 28 21:39:13 2007 INFO Continue...: # 14:06:15.0 - 14:16:25.0 21 3 3C129 [0, 1] # 14:20:45.0 - 14:22:25.0 22 2 0420+417 [0, 1] # 14:45:15.0 - 14:55:25.0 23 3 3C129 [0, 1] # 14:59:45.0 - 15:01:25.0 24 2 0420+417 [0, 1] # 15:18:25.0 - 15:25:35.0 25 3 3C129 [0, 1] # 15:30:05.0 - 15:31:45.0 26 2 0420+417 [0, 1] # 15:55:05.0 - 16:05:15.0 27 3 3C129 [0, 1] # 16:10:05.0 - 16:11:45.0 28 2 0420+417 [0, 1] # 16:34:55.0 - 16:45:05.0 29 3 3C129 [0, 1] # 16:49:35.0 - 16:51:15.0 30 2 0420+417 [0, 1] # 17:06:35.0 - 17:09:45.0 31 4 0518+165 [0, 1] # 17:21:45.0 - 17:24:55.0 32 5 0134+329 [0, 1] # # Wed Nov 28 21:39:13 2007 INFO listobs::ms::summary: # ObservationID = 2 ArrayID = 0 # Date Timerange Scan FldId FieldName SpwIds # 03-Nov-1994/04:05:15.0 - 04:07:25.0 33 5 0134+329 [0, 1] # 07:33:35.0 - 07:35:15.0 34 2 0420+417 [0, 1] # 07:42:35.0 - 07:45:45.0 35 4 0518+165 [0, 1] # 08:18:05.0 - 08:23:25.0 36 3 3C129 [0, 1] # 08:27:35.0 - 08:51:05.0 37 2 0420+417 [0, 1] # 09:07:05.0 - 09:12:25.0 38 3 3C129 [0, 1] # 09:17:05.0 - 09:18:45.0 39 2 0420+417 [0, 1] # 09:40:25.0 - 09:45:05.0 40 3 3C129 [0, 1] # 09:49:35.0 - 09:51:45.0 41 2 0420+417 [0, 1] # # Wed Nov 28 21:39:14 2007 INFO listobs::ms::summary: # Fields: 6 # ID Code Name Right Ascension Declination Epoch # 0 A 2345-167 23:45:27.68 -16.47.52.60 B1950_VLA # 1 A2597 23:22:43.70 -12.23.56.00 B1950_VLA # 2 B 0420+417 04:20:27.94 +41.43.08.04 B1950_VLA # 3 3C129 04:45:31.69 +44.55.19.95 B1950_VLA # 4 C 0518+165 05:18:16.53 +16.35.26.90 B1950_VLA # 5 C 0134+329 01:34:49.83 +32.54.20.52 B1950_VLA # # Wed Nov 28 21:39:14 2007 INFO listobs::ms::summary: # Spectral Windows: (2 unique spectral windows and 1 unique polarization setups) # SpwID #Chans Frame Ch1(MHz) ChanWid(kHz)TotBW(kHz) Ref(MHz) Corrs # 0 1 TOPO 4885.1 50000 50000 4885.1 RR RL LR LL # 1 1 TOPO 4585.1 50000 50000 4585.1 RR RL LR LL # # Wed Nov 28 21:39:14 2007 INFO listobs::ms::summary: # Feeds: 47: printing first row only # Antenna Spectral Window # Receptors Polarizations # 1 -1 2 [ R, L] # # Wed Nov 28 21:39:14 2007 INFO listobs::ms::summary: # Antennas: 47: # ID Name Station Diam. Long. Lat. # 0 VA02 VLA:W20 25.0 m -107.38.21.4 +33.53.19.5 # 1 VA28 VLA:W16 25.0 m -107.37.57.4 +33.53.33.0 # 2 VA08 VLA:W8 25.0 m -107.37.21.6 +33.53.53.0 # 3 VA03 VLA:W12 25.0 m -107.37.37.4 +33.53.44.2 # 4 VA12 VLA:W4 25.0 m -107.37.10.8 +33.53.59.1 # 5 VA18 VLA:W32 25.0 m -107.39.54.8 +33.52.27.2 # 6 VA10 VLA:W28 25.0 m -107.39.20.2 +33.52.46.6 # 7 VA20 VLA:W36 25.0 m -107.40.32.6 +33.52.06.0 # 8 VA21 VLA:W24 25.0 m -107.38.49.0 +33.53.04.0 # 9 VA17 VLA:E20 25.0 m -107.35.43.6 +33.53.29.9 # 10 VA07 VLA:E16 25.0 m -107.36.09.8 +33.53.40.0 # 11 VA04 VLA:E8 25.0 m -107.36.48.9 +33.53.55.1 # 12 VA16 VLA:E12 25.0 m -107.36.31.7 +33.53.48.5 # 13 VA22 VLA:E4 25.0 m -107.37.00.8 +33.53.59.7 # 14 VA23 VLA:E32 25.0 m -107.34.01.5 +33.52.50.3 # 15 VA05 VLA:E28 25.0 m -107.34.39.3 +33.53.04.9 # 16 VA06 VLA:E36 25.0 m -107.33.20.2 +33.52.34.3 # 17 VA24 VLA:E24 25.0 m -107.35.13.4 +33.53.18.1 # 18 VA11 VLA:N24 25.0 m -107.37.16.1 +33.55.37.7 # # Wed Nov 28 21:39:14 2007 INFO Continue...: # 19 VA15 VLA:N20 25.0 m -107.37.13.2 +33.55.09.5 # 20 VA27 VLA:N8 25.0 m -107.37.07.5 +33.54.15.8 # 21 VA13 VLA:N12 25.0 m -107.37.09.0 +33.54.30.0 # 22 VA25 VLA:N16 25.0 m -107.37.10.9 +33.54.48.0 # 23 VA14 VLA:N4 25.0 m -107.37.06.5 +33.54.06.1 # 24 VA26 VLA:N36 25.0 m -107.37.25.6 +33.57.07.6 # 25 VA01 VLA:N32 25.0 m -107.37.22.0 +33.56.33.6 # 26 VA19 VLA:N28 25.0 m -107.37.18.7 +33.56.02.5 # 27 VA02 VLA:W10 25.0 m -107.37.28.9 +33.53.48.9 # 28 VA03 VLA:W2 25.0 m -107.37.07.4 +33.54.00.9 # 29 VA18 VLA:W14 25.0 m -107.37.46.9 +33.53.38.9 # 30 VA10 VLA:W12 25.0 m -107.37.37.4 +33.53.44.2 # 31 VA20 VLA:W18 25.0 m -107.38.08.9 +33.53.26.5 # 32 VA21 VLA:W6 25.0 m -107.37.15.6 +33.53.56.4 # 33 VA17 VLA:E8 25.0 m -107.36.48.9 +33.53.55.1 # 34 VA04 VLA:E2 25.0 m -107.37.04.4 +33.54.01.1 # 35 VA16 VLA:E6 25.0 m -107.36.55.6 +33.53.57.7 # 36 VA23 VLA:E14 25.0 m -107.36.21.3 +33.53.44.5 # 37 VA24 VLA:E10 25.0 m -107.36.40.9 +33.53.52.0 # 38 VA06 VLA:E18 25.0 m -107.35.57.2 +33.53.35.1 # 39 VA05 VLA:E12 25.0 m -107.36.31.7 +33.53.48.5 # # Wed Nov 28 21:39:14 2007 INFO Continue...: # 40 VA11 VLA:N12 25.0 m -107.37.09.0 +33.54.30.0 # 41 VA15 VLA:N10 25.0 m -107.37.08.2 +33.54.22.4 # 42 VA13 VLA:N2 25.0 m -107.37.06.2 +33.54.03.5 # 43 VA25 VLA:N6 25.0 m -107.37.06.9 +33.54.10.3 # 44 VA26 VLA:N18 25.0 m -107.37.12.0 +33.54.58.3 # 45 VA01 VLA:N16 25.0 m -107.37.10.9 +33.54.48.0 # 46 VA09 VLA:N14 25.0 m -107.37.09.9 +33.54.38.5 # # Wed Nov 28 21:39:14 2007 INFO listobs::ms::summary: # Tables: # MAIN 677430 rows # ANTENNA 47 rows # DATA_DESCRIPTION 2 rows # DOPPLER 2 rows # FEED 47 rows # FIELD 6 rows # FLAG_CMD # FREQ_OFFSET # HISTORY 6 rows # OBSERVATION 3 rows # POINTING # POLARIZATION 1 row # PROCESSOR # SOURCE 6 rows # SPECTRAL_WINDOW 2 rows # STATE # SYSCAL # WEATHER # # # #===================================================================== # Data Examination and Flagging #===================================================================== # # Get rid of the autocorrelations from the MS # print '--Flagautocorr--' print "You can use flagautocorr to zap auto-correlations if necessary." print "By default, importvla does not load them, so we don't need to run it." # Don't default this one either # flagautocorr() # #===================================================================== # # Use Plotxy to interactively flag the data # print '--Plotxy--' default('plotxy') print "Now we use plotxy to examine and interactively flag data." print "We'll just look at parallel-hand correlations here, but you can look at the" print "cross-hand polarizations by setting correlation=\"RL LR\"" vis = msfile # The fields we are interested in: 0134+329 (3C48), 0518+165 (3C138), 0420+417 (phase cal), 3C129 selectdata = True # Plot only the RR and LL for now correlation = 'RR LL' # Plot amplitude vs. uvdist xaxis = 'uvdist' yaxis = 'amp' multicolor = 'both' # The easiest thing is to iterate over antennas iteration = 'antenna' print "" print "Showing 0134+329 with iteration='antenna' " print "Use Next button to step through antennas" print "These data look pretty good to me... I wouldn't do anything." # Do the primary calibrators first field = '0134+329' plotxy() # Pause script if you are running in scriptmode if scriptmode: user_check=raw_input('Return to continue script\n') print "" print "Showing 0518+165 with iteration='antenna' " print "Use Next button to step through antennas" print "These data look pretty good to me... I wouldn't do anything." field = '0518+165' plotxy() # Pause script if you are running in scriptmode if scriptmode: user_check=raw_input('Return to continue script\n') # Now do the phase calibrator print "" print "Plotting 0420+417 RR LL all antennas" print "This is good data, but if you feel the need to flag something..." print "Mark a box around a bit of it and hit Locate" print "Look in logger to see what it is" # Pause script if you are running in scriptmode if scriptmode: user_check=raw_input('Return to continue script\n') plotxy() # Pause script if you are running in scriptmode if scriptmode: user_check=raw_input('Return to continue script\n') #--------------------------------------------------------------------- # Finally, do 3C129 field = '3C129' correlation = 'RR LL' iteration = '' xaxis = 'uvdist' plotxy() # Pause script if you are running in scriptmode if scriptmode: user_check=raw_input('Return to continue script\n') print "" print "Now plotting vs. time, one antenna at a time" xaxis = 'time' iteration = 'antenna' plotxy() print "" print "Step through antennas with Next" print "" print "You can draw little boxes around the outliers and Flag" print "You can also use Locate to find where they come from" print "" print "You will be drawing many tiny boxes, so remember you can" print "use the ESC key to get rid of the most recent box if you" print "make a mistake." # Pause script if you are running in scriptmode if scriptmode: user_check=raw_input('Return to continue script\n') print "Done with plotxy!" # #===================================================================== # # Use Flagmanager to save a copy of the flags so far # print '--Flagmanager--' default('flagmanager') print "Now will use flagmanager to save a copy of the flags we just made" print "These are named basic_flags" vis = msfile mode = 'save' versionname = 'basic flags' comment = 'Plotxy flags' merge = 'replace' flagmanager() # #===================================================================== # # You can use Flagdata to explicitly clip the data also # #print '--Flagdata--' #default('flagdata') # #print "As a demonstration show how to clip the data with flagdata" #print "Note we had already flagged many of these interactively" # #vis = msfile # ## Set some clipping regions #mode = 'manualflag' #clipcolumn = 'DATA' #clipoutside = False # ## Clip calibraters #field = '1331+305' #clipexpr = 'ABS RR' #clipminmax = [0.0,0.75] #flagdata() # #clipexpr = 'ABS LL' #clipminmax = [0.0,0.75] #flagdata() # #clipexpr = 'ABS RL' #clipminmax = [0.0,0.055] #flagdata() # #clipexpr = 'ABS LR' #clipminmax = [0.0,0.055] #flagdata() # #field = '0137+331' #clipexpr = 'ABS RR' #clipminmax = [0.0,0.55] #flagdata() # #clipexpr = 'ABS LL' #clipminmax = [0.0,0.55] #flagdata() # ## You can also do the antenna edits on 0137+331 and JUPITER ## with flagdata # ## Done with flagging # #===================================================================== # # Use Flagmanager to list all saved versions # print '--Flagmanager--' default('flagmanager') print "Now will use flagmanager to list the versions we saved" vis = msfile mode = 'list' flagmanager() # #===================================================================== # Calibration #===================================================================== # # Set the fluxes of the primary calibrator(s) # print '--Setjy--' default('setjy') print "Use setjy to set flux of 0134+329 (3C48) and 0518+165 (3C138)" vis = msfile # # 1331+305 = 3C286 is our primary calibrator # Setjy knows about these sources so we only need to set the field field = '0134+329' setjy() field = '0518+165' setjy() # # You should see something like this in the logger and casapy.log file: # # 1331+305 spwid= 0 [I=7.462, Q=0, U=0, V=0] Jy, (Perley-Taylor 99) # 1331+305 spwid= 1 [I=7.51, Q=0, U=0, V=0] Jy, (Perley-Taylor 99) # print "Look in logger for the fluxes" # #===================================================================== # # Initial gain calibration # print '--Gaincal--' default('gaincal') print "Solve for antenna gains on the calibrators:" print "0134+329 (3C48), 0518+165 (3C138), and 0420+417" print "We have 2 single-channel continuum spw" print "Do not want bandpass calibration" vis = msfile # set the name for the output gain caltable gtable = prefix + '.gcal' caltable = gtable print "Output gain cal table will be "+gtable # We have 2 IFs (SPW 0,1) with one channel each # selection is via the field and spw strings field = '0134+329,0518+165,0420+417' spw = '' # a-priori calibration application # atmospheric optical depth (turn off) gaincurve = True opacity = 0.0 # scan-based G solutions for both amplitude and phase gaintype = 'G' solint = 0. calmode = 'ap' # reference antenna 11 (11=VLA:N1) refant = '11' # minimum SNR 3 minsnr = 3 gaincal() # #===================================================================== # # Bootstrap flux scale # print '--Fluxscale--' default('fluxscale') print "Use fluxscale to rescale the gains. We'll do this in-place." vis = msfile # point to our first gain cal table caltable = gtable # we will be using the sources we ran setjy on as # our flux standard references reference = '0134+329,0518+165' # we want to transfer the flux to our other gain cal source 0137+331 # to bring its gain amplitues in line with the absolute scale transfer = '0420+417' fluxscale() # You should see in the logger something like: #Flux density for 0137+331 in SpW=0 is: # 5.42575 +/- 0.00285011 (SNR = 1903.7, nAnt= 27) #Flux density for 0137+331 in SpW=1 is: # 5.46569 +/- 0.00301326 (SNR = 1813.88, nAnt= 27) #===================================================================== # # Interpolate the gains onto all targets # print '--Accum--' default('accum') print "This will interpolate the gains onto 3C129" vis = msfile tablein = '' incrtable = gtable calfield = '0134+329,0518+165,0420+417' # set the name for the output interpolated caltable atable = prefix + '.accum' caltable = atable print "Output cumulative gain table will be "+atable # linear interpolation interp = 'linear' # make 10s entries... accumtime = 10.0 accum() #===================================================================== # # Correct the data # (This will put calibrated data into the CORRECTED_DATA column) # print '--ApplyCal--' default('applycal') print "This will apply the calibration to the DATA" print "Fills CORRECTED_DATA" vis = msfile # Start with the interpolated fluxscale/gain table gaintable = atable # Since we did gaincurve=True in gaincal, we need it here also gaincurve = True opacity=0.0 # select all fields field = '' spw = '' selectdata = False # do not need to select subset since we did accum # (note that correct only does 'nearest' interp) gainselect = '' applycal() # #===================================================================== # # Now split the target data # print '--Split--' default('split') vis = msfile # Now we write out the corrected data to a new MS # Make an output vis file srcsplitms = prefix + '.split.ms' outputvis = srcsplitms print "Split 3C129 data into new ms "+srcsplitms # Select the field field = '3C129' spw = '' # pick off the CORRECTED_DATA column datacolumn = 'corrected' split() #===================================================================== # # Export the 3C129 data as UVFITS # Start with the split file. # print '--Export UVFITS--' default('exportuvfits') srcuvfits = prefix + '.split.uvfits' print "Writing split data to UVFITS file "+srcuvfits vis = srcsplitms fitsfile = srcuvfits # Since this is a split dataset, the calibrated data is # in the DATA column already. datacolumn = 'data' # Write as a multisource UVFITS (with SU table) # even though it will have only one field in it multisource = True # Run asynchronously so as not to interfere with other tasks # (BETA: also avoids crash on next importuvfits) async = True exportuvfits() # #===================================================================== # FIRST CLEAN / SELFCAL CYCLE #===================================================================== # # Now clean an image of 3C129 # print '--Clean 1--' default('clean') # Pick up our split source data vis = srcsplitms # Make an image root file name imname1 = prefix + '.clean1' imagename = imname1 print "Output images will be prefixed with "+imname1 # Set up the output continuum image (single plane mfs) mode = 'mfs' stokes = 'I' print "Will be a single MFS continuum image" # NOTE: current version field='' doesnt work field = '3C129' # Combine all spw spw = '' # This is B and C-config VLA 6cm (4.85GHz) obs # Check the observational status summary # Primary beam FWHM = 45'/f_GHz = 557" # Synthesized beam FWHM = 1" # RMS in 10min (600s) = 0.06 mJy (thats now, but close enough) # Set the output image size and cell size (arcsec) # .3" will give 3x oversampling # 3700 pix will cover 2x Primary Beam # we'll use 2048 alg = 'clark' imsize = [2048,2048] cell = ['0.3arcsec','0.3arcsec'] # NOTE: will eventually have an imadvise task to give you this # information # Standard gain factor 0.1 gain = 0.1 # Fix maximum number of iterations niter = 100000 # Note - we can change niter and threshold interactively # during clean # Set up the weighting # Use Briggs weighting (a moderate value, on the uniform side) weighting = 'briggs' rmode = 'norm' robust = 0.5 # No clean mask mask = '' # Use interactive clean mode cleanbox = 'interactive' # Moderate number of iter per interactive cycle npercycle = 100 clean() # When the interactive clean window comes up, use the right-mouse # to draw rectangles around obvious emission double-right-clicking # inside them to add to the flag region. You can also assign the # right-mouse to polygon region drawing by right-clicking on the # polygon drawing icon in the toolbar. When you are happy with # the region, click 'Done Flagging' and it will go and clean another # 100 iterations. When done, click 'Stop'. # Set up variables clnimage1 = imname1+'.image' clnmodel1 = imname1+'.model' clnresid1 = imname1+'.residual' clnmask1 = imname1+'.clean_interactive.mask' print "" print "Final clean model is "+clnmodel1 print "Final restored clean image is "+clnimage1 print "The clean residual image is "+clnresid1 print "Your final clean mask is "+clnmask1 print "" print "This is the final restored clean image in the viewer" print "Zoom in and set levels to see faint emission" print "Use rectangle drawing tool to box off source" print "Double-click inside to print statistics" print "Move box on-source and get the max" print "Calcualte DynRange = MAXon/RMSoff" print "I got 0.036 / 9.7E-5 = 370 print "Still not as good as it can be - lets selfcal" print "Close viewer panel when done" # #--------------------------------------------------------------------- # # If you did not do interactive clean, bring up viewer manually viewer(clnimage1,'image') # Pause script if you are running in scriptmode if scriptmode: user_check=raw_input('Return to continue script\n') # You can use the right-mouse to draw a box in the lower right # corner of the image away from emission, the double-click inside # to bring up statistics. Use the right-mouse to grab this box # and move it up over the source and double-click again. Look at # the terminal output... # #--------------------------------------------------------------------- # # Self-cal using clean model # # Note: clean will have left FT of model in the MODEL_DATA column # If you've done something in between, can use the ft task to # do this manually. # print '--SelfCal 1--' default('gaincal') vis = srcsplitms print "Will self-cal using MODEL_DATA left in MS by clean" # New gain table selfcaltab1 = srcsplitms + '.selfcal1' caltable = selfcaltab1 print "Will write gain table "+selfcaltab1 # Don't need a-priori cals selectdata = False gaincurve = False opacity = 0.0 # This choice seemed to work refant = '11' # Lets do phase-only first time around gaintype = 'G' calmode = 'p' # Do scan-based solutions with SNR>3 solint = 0.0 minsnr = 3.0 # Do not need to normalize (let gains float) solnorm = False gaincal() # #--------------------------------------------------------------------- # # Correct the data (no need for interpolation this stage) # print '--ApplyCal--' default('applycal') vis = srcsplitms print "Will apply self-cal table to over-write CORRECTED_DATA in MS" gaintable = selfcaltab1 gaincurve = False opacity = 0.0 field = '' spw = '' selectdata = False calwt = True applycal() # Self-cal is now in CORRECTED_DATA column of split ms # #===================================================================== # SECOND CLEAN / SELFCAL CYCLE #===================================================================== # print '--Clean 2--' default('clean') print "Now clean on self-calibrated data" vis = srcsplitms imname2 = prefix + '.clean2' imagename = imname2 field = '*' spw = '' mode = 'mfs' gain = 0.1 niter = 100000 alg = 'clark' imsize = [2048,2048] cell = ['0.3arcsec','0.3arcsec'] weighting = 'briggs' rmode = 'norm' robust = 0.5 cleanbox = 'interactive' npercycle = 100 clean() # Set up variables clnimage2 = imname2+'.image' clnmodel2 = imname2+'.model' clnresid2 = imname2+'.residual' clnmask2 = imname2+'.clean_interactive.mask' print "" print "Final clean model is "+clnmodel2 print "Final restored clean image is "+clnimage2 print "The clean residual image is "+clnresid2 print "Your final clean mask is "+clnmask2 print "" print "This is the final restored clean image in the viewer" print "Zoom in and set levels to see faint emission" print "Use rectangle drawing tool to box off source" print "Double-click inside to print statistics" print "Move box on-source and get the max" print "Calcualte DynRange = MAXon/RMSoff" print "This time I got 0.036 / 9.1E-5 = 395 (not much better)" print "Still not as good as it can be - lets selfcal again" print "Close viewer panel when done" # #--------------------------------------------------------------------- # # If you did not do interactive clean, bring up viewer manually viewer(clnimage1,'image') # Pause script if you are running in scriptmode if scriptmode: user_check=raw_input('Return to continue script\n') # Note that the dynamic range will depend on how deep you take # the interactive clean and how you draw the box for the stats. # #--------------------------------------------------------------------- # # Next self-cal cycle # print '--SelfCal 2--' default('gaincal') print "Self-cal again using MODEL_DATA from second clean" vis = srcsplitms selfcaltab2 = srcsplitms + '.selfcal2' caltable = selfcaltab2 print "Write self-cal table "+selfcaltab2 selectdata = False gaincurve = False opacity = 0.0 refant = '11' # This time amp+phase on 10s timescales SNR>1 gaintype = 'G' calmode = 'ap' solint = 10.0 minsnr = 1.0 solnorm = False gaincal() # # It is useful to put this up in plotcal # #--------------------------------------------------------------------- # print '--PlotCal--' default('plotcal') caltable = selfcaltab2 multiplot = True yaxis = 'amp' plotcal() # Use the Next button to iterate over antennas print "" print "Looking at amplitude in self-cal table "+caltable print "Note coherence of solutions" # Pause script if you are running in scriptmode if scriptmode: user_check=raw_input('Return to continue script\n') yaxis = 'phase' plotcal() # # You can see it is not too noisy. print "" print "Now look at phases" # # Pause script if you are running in scriptmode if scriptmode: user_check=raw_input('Return to continue script\n') # Lets do some smoothing anyway. # #--------------------------------------------------------------------- # # Smooth calibration solutions # print '--Smooth--' default('smoothcal') print "Lets do smoothcal anyway as a demonstration" vis = srcsplitms tablein = selfcaltab2 smoothcaltab2 = srcsplitms + '.smoothcal2' caltable = smoothcaltab2 print "Will write smoothed table "+smoothcaltab2 # Do a 30s boxcar average smoothtype = 'mean' smoothtime = 30.0 smoothcal() # If you put into plotcal you'll see the results # For example, you can grap the inputs from the last # time you ran plotcal, set the new tablename, and plot! #run plotcal.last #tablein = smoothcaltab2 #plotcal() # #--------------------------------------------------------------------- # # Correct the data # print '--ApplyCal--' default('applycal') print "Calibrate data using the smoothed table" vis = srcsplitms gaintable = smoothcaltab2 gaincurve = False opacity = 0.0 field = '' spw = '' selectdata = False applycal() # #===================================================================== # THIRD CLEAN / SELFCAL CYCLE #===================================================================== # print '--Clean 3--' default('clean') print "This is the 3rd clean using smoothed 2nd self-cal" vis = srcsplitms imname3 = prefix + '.clean3' imagename = imname3 field = '*' spw = '' mode = 'mfs' gain = 0.1 niter = 100000 alg = 'clark' imsize = [2048,2048] cell = ['0.3arcsec','0.3arcsec'] weighting = 'briggs' rmode = 'norm' robust = 0.5 cleanbox = 'interactive' npercycle = 100 clean() # Cleans alot deeper # You can change the npercycle to larger numbers # (like 250 or so) as you get deeper also. # Set up variables clnimage3 = imname3+'.image' clnmodel3 = imname3+'.model' clnresid3 = imname3+'.residual' clnmask3 = imname3+'.clean_interactive.mask' print "" print "Final clean model is "+clnmodel3 print "Final restored clean image is "+clnimage3 print "The clean residual image is "+clnresid3 print "Your final clean mask is "+clnmask3 print "" print "This is the final restored clean image in the viewer" print "Zoom in and set levels to see faint emission" print "Use rectangle drawing tool to box off source" print "Double-click inside to print statistics" print "Move box on-source and get the max" print "Calcualte DynRange = MAXon/RMSoff" print "" print "Close viewer panel when done" # #--------------------------------------------------------------------- # # If you did not do interactive clean, bring up viewer manually viewer(clnimage1,'image') # Pause script if you are running in scriptmode if scriptmode: user_check=raw_input('Return to continue script\n') # Note that the dynamic range you get will depend on how deep you # take the interactive clean and how you draw the box for the stats. # # This will probably take several more careful self-cal cycles. # Set up final variables clnimage = clnimage3 clnmodel = clnmodel3 clnresid = clnresid3 clnmask = clnmask3 print "" print "After this script is done you can continue on with" print "more self-cal, or try different cleaning options" #===================================================================== # # Export the Final CLEAN Image as FITS # print '--Final Export CLEAN FITS--' default('exportfits') clnfits = prefix + '.clean.fits' print "Write out final FITS restored clean image "+clnfits imagename = clnimage fitsimage = clnfits # Run asynchronously so as not to interfere with other tasks # (BETA: also avoids crash on next importfits) async = True exportfits() #===================================================================== # # Export the Final Self-Calibrated data as UVFITS # print '--Final Export UVFITS--' default('exportuvfits') caluvfits = prefix + '.selfcal.uvfits' print "Write out final UVFITS self-cal data MS "+caluvfits vis = srcsplitms fitsfile = caluvfits # The self-calibrated data is in the CORRECTED_DATA column datacolumn = 'corrected' # Write as a multisource UVFITS (with SU table) # even though it will have only one field in it multisource = True # Run asynchronously so as not to interfere with other tasks # (BETA: also avoids crash on next importuvfits) async = True exportuvfits() # #===================================================================== # Image Analysis #===================================================================== # # Can do some image statistics if you wish # Treat this like a regression script ## WARNING: currently requires toolkit ## #print "Using MS to get stats of UV data" #print "Using ia() tool to get stats of clean image" #print "" # #print ' Results ' #print ' =============== ' # #print '' ## Pull the max src amp value out of the MS #ms.open(srcsplitms) #thistest_src = max(ms.range(["amplitude"]).get('amplitude')) #oldtest_src = 4.92000198364 #print ' MS max amplitude should be ',oldtest_src #print ' Found : Max in MS = ',thistest_src #diff_src = abs((oldtest_src-thistest_src)/oldtest_src) #print ' Difference (fractional) = ',diff_src # #ms.close() # #print '' ## Pull the max and rms from the clean image #ia.open(clnimage) #on_statistics=ia.statistics() #thistest_immax=on_statistics['max'][0] #oldtest_immax = 1.07732224464 #print ' Clean image ON-SRC max should be ',oldtest_immax #print ' Found : Max in image = ',thistest_immax #diff_immax = abs((oldtest_immax-thistest_immax)/oldtest_immax) #print ' Difference (fractional) = ',diff_immax # #print '' ## Now do stats in the lower right corner of the image #box = ia.setboxregion([0.75,0.00],[1.00,0.25],frac=true) #off_statistics=ia.statistics(region=box) #thistest_imrms=off_statistics['rms'][0] #oldtest_imrms = 0.0010449 #print ' Clean image OFF-SRC rms should be ',oldtest_imrms #print ' Found : rms in image = ',thistest_imrms #diff_imrms = abs((oldtest_imrms-thistest_imrms)/oldtest_imrms) #print ' Difference (fractional) = ',diff_imrms # #print '' #print ' Final Clean image Dynamic Range = ',thistest_immax/thistest_imrms #print '' #print ' =============== ' # #ia.close() # #print '' #print '--- Done ---' # ## ##=====================================================================