import os from recipes.almapolhelpers import * # BeamCrossRedux.py # # Written by Chat Hull, Jul 2015 # Final updates for data release: Apr 2020 ##### HOW TO RUN THE CODE IN CASA (VERSION <= 5.6.1) # (1) Place this script in a directory, and create a subdirectory for the data # (2) Copy raw ASDMs into the data subdirectory # (3) cd into the data subdirectory # Then: # execfile('../beamCrossRedux.py') # setup(B3_3) # S = onAxisRedux(B3_3) # offAxisRedux_onAxisCal(B3_3, S) # offAxisImageAnalysis(B3_3) # Above, replace B3_3 with whichever trial you want to reduce # Details for different 11x11 trials # B3-3: 31 ants, refant = DV21 -- 06mar15 # B5-1: 9 ants, refant = DV21 -- 10nov16 # B6-1: 32 ants, refant = DA59 -- 08may15 # B7-1: 35 ants, refant = DA59 -- 09may15 B3_3 = { 'source' : '3c279', 'refant' : 'DV21', 'datasets' : ['uid___A002_X9b98ec_X94d.ms', 'uid___A002_X9b98ec_Xb98.ms', 'uid___A002_X9b98ec_Xde3.ms'], 'flags' : [], 'onAxisScans' : '100,161,222', 'KcrsScans' : '222,223,284,345,346', 'width' : [64,64,64,64], 'size' : '11x11', 'imsize' : [300], 'band' : '3', 'dir' : 'B3-3', } B5_1 = { 'source' : '3c279', 'refant' : 'DV21', 'datasets' : ['uid___A002_Xba6edb_X31db.ms', 'uid___A002_Xba6edb_X3988.ms', 'uid___A002_Xba6edb_X3d57.ms'], 'flags' : [], 'onAxisScans' : '100,161,222', 'KcrsScans' : '284,345,346,407,468', 'width' : [64,64,64,64], 'size' : '11x11', 'imsize' : [300], 'band' : '5', 'dir' : 'B5-1', } B6_1 = { 'source' : '3c279', 'refant' : 'DA59', 'datasets' : ['uid___A002_Xa018c4_X743.ms', 'uid___A002_Xa018c4_Xa3d.ms', 'uid___A002_Xa018c4_Xdc7.ms'], 'flags' : ['flagdata(vis=vis, antenna="PM03")', 'flagdata(vis=vis, antenna="PM04,DA45", spw="3")', 'flagdata(vis=vis, antenna="PM04", spw="3")'], 'onAxisScans' : '100,161,222', 'KcrsScans' : '100,161,222,223,284', 'width' : [64,64,64,64], 'size' : '11x11', 'imsize' : [300], 'band' : '6', 'dir' : 'B6-1', } B7_1 = { 'source' : '3c279', 'refant' : 'DA59', 'datasets' : ['uid___A002_Xa018c4_X3b92.ms', 'uid___A002_Xa018c4_X3de3.ms', 'uid___A002_Xa018c4_X402a.ms'], 'flags' : ['flagdata(vis=vis, antenna="PM03")', 'flagdata(vis=vis, antenna="PM04,DA45", spw="3")', 'flagdata(vis=vis, antenna="PM04", spw="3")'], 'onAxisScans' : '100,161,222', 'KcrsScans' : '100,161,222,223,284,345,346', 'width' : [64,64,64,64], 'size' : '11x11', 'imsize' : [300], 'band' : '7', 'dir' : 'B7-1', } def setup(trial): # Check directory if os.path.abspath('.').split('/')[-1] != trial['dir'] : raise Exception("You're not in the correct directory!") # Get refant refant = trial['refant'] # Get datasets datasets = trial['datasets'] # Import ASDMs for vis in trial['datasets']: os.system('rm -rf '+vis) os.system('rm -rf '+vis+'.flagversions') importasdm(asdm=vis.strip('.ms'), vis=vis) # Solve degeneracy of subscans, and assign a new scan number to each field # This splits the 1 original scan into its 222 individual scans for vis in trial['datasets']: tb.open(vis,nomodify=False) scno=tb.getcol('SCAN_NUMBER') scno*=100 # makes nominal scan numbers far apart stid=tb.getcol('STATE_ID') scno+=stid # add state_id to scan_number tb.putcol('SCAN_NUMBER',scno) tb.close() #### FLAGGING for vis in trial['datasets'] : #### General flagging for all datasets # Flag edge channels flagdata(vis=vis, spw='*FULL_RES*:0~4;59~63') # Flag shadowed antennas flagdata(vis=vis, mode='shadow', flagbackup=False) # Flag autocorrelations flagdata(vis=vis, mode='manual', autocorr=True, flagbackup=False) #### Specific flagging for each dataset for f in trial['flags'] : exec(f) #### CONCAT AND SPLIT DATA # Concatenate 3 SBs os.system('rm -rf concat.ms') concat(vis=datasets,concatvis='concat.ms') # Channel-average concatenated data and split out appropriate spws os.system('rm -rf concat.ms.split') os.system('rm -rf concat.ms.split.flagversions') split2(vis='concat.ms', outputvis='concat.ms.split', spw='*FULL_RES*', datacolumn='data') # Make backup of full dataset os.system('cp -r concat.ms.split concat.ms.split.backup') #### ON-AXIS PREP # Split out on-axis field (first/middle/last scans) from each of the 3 datasets for vis in trial['datasets'] : os.system('rm -rf {0}.onaxis.split'.format(vis)) split2(vis=vis, outputvis='{0}.onaxis.split'.format(vis), scan=trial['onAxisScans'], spw='*FULL_RES*', datacolumn='data') # Concat 3 on-axis observations onaxis_vis=[datasets[0]+'.onaxis.split', datasets[1]+'.onaxis.split', datasets[2]+'.onaxis.split'] os.system('rm -rf onaxis_raw.ms') concat(vis=onaxis_vis, concatvis='onaxis_raw.ms') # listobs os.system('rm -rf onaxis.ms.listobs') listobs(vis='onaxis_raw.ms',listfile='onaxis.ms.listobs') os.system('rm -rf concat.ms.split.listobs') listobs(vis='concat.ms.split',listfile='concat.ms.split.listobs') def onAxisRedux(trial) : # Check directory if os.path.abspath('.').split('/')[-1] != trial['dir'] : raise Exception("You're not in the correct directory!") # Get refant refant = trial['refant'] # Get datasets datasets = trial['datasets'] # Decide whether to pop plots or not interact = False # Copy original data file os.system('rm -rf onaxis.ms') os.system('cp -r onaxis_raw.ms onaxis.ms') # On-axis data msdata='onaxis.ms' msdata_cal='onaxis.cal.ms' # We would normally do WVR and tsys correction at this point, but the Tsys weren't in all datasets ########################### ## Bandpass calibration ### ########################### # In preparation to the bandpass calibration, we correct the phase time variation. # Do not choose edge channels and atmospheric absorption # This would normally be taken care of in the Tsys data, but we don't have Tsys # This is the initial, pre-bandpass phase calibration os.system('rm -rf '+msdata+'.G0ph') gaincal(vis=msdata, caltable=msdata+'.G0ph', spw='*', gaintype='G', solint='int', # Solve for each integration so that we can combine in next step calmode='p', refant=refant, smodel=[1,0,0,0], # This is the default interp='nearest') # Bandpass calibration using G0ph os.system('rm -rf '+msdata+'.Bpass') bandpass(vis=msdata, caltable=msdata+'.Bpass', solint='inf', # combines all integrations combine='scan', # Gives 3 separate bandpass solutions, because of 3 separate obs blocks refant=refant, solnorm=True, # Normalizes each bandpass solution to 1 gaintable=[msdata+'.G0ph'], interp=['nearest']) # Plot bandpass table if interact : os.system('rm -rf ./plotbandpass/') os.system('mkdir ./plotbandpass/') plotbandpass(caltable=msdata+'.Bpass', xaxis='chan', yaxis='both', figfile='./plotbandpass/BpassPlot.png', figfileSequential=True, interactive=False) ################################ # Gain (amp&phase) calibration # ################################ os.system('rm -rf '+msdata+'.G1') gaincal(vis=msdata, caltable=msdata+'.G1', # default gaintype=G, default calmode='ap' solint='int', # one for each integration refant=refant, smodel=[1,0,0,0], gaintable=[msdata+'.Bpass'], # Could keep using .G0, but don't here; .G0 was used to get good .Bpass interp=['nearest']) ################################################### # Applying bandpass and gain cal (just for check) # ################################################### applycal(vis=msdata, calwt=False, # =True only for Tsys gaintable=[msdata+'.Bpass',msdata+'.G1'], interp=['nearest','linear']) # Use 'nearest' when fields are different (cal vs. source) # Use 'linear' when you're calibrating the source itself with its own observations ################################# # QU estimate from G1 # ################################# qu=qufromgain(msdata+'.G1') # This is because we have crossed-linear feeds! # We assumed unpolarized source; thus, gains compensate for actual polarization # EXAMPLE TEXT: # Latitude = -23.0290058252 # Found as many as 1 fields. # Found as many as 1 spws. # Fld= 0 Spw= 0 (B=07, PA offset=-52.5deg) Gx/Gy= 1.36479084722 Q= 0.160990495326 U= 0.0704752589156 P= 0.175740438444 X= 11.8209489956 # For field id = 0 there are 1 good spws. # Spw mean: Fld= 0 Q= 0.160990495326 U= 0.0704752589156 (rms= 0.0 0.0 ) P= 0.175740438444 X= 11.8209489956 ################################### ## Cross-hand delay calibration ### ################################### # Effectively a bandpass on the XY and YX # NEED TO CHANGE SCAN NUMBERS FOR DIFFERENT DATASETS # KCROSS scans # Use the scans where the polarization ratio is changing fastest # Chose the middle three scans (first and third are pairs of scans) # To choose scan, look for scan that is on the max slope of poln-ratio-gain vs. scan if interact: plotcal('onaxis.ms.G1', 'scan', 'amp', field='', poln='/') os.system('rm -rf '+msdata+'.Kcrs') gaincal(vis=msdata, caltable=msdata+'.Kcrs', selectdata=True, # This is so we can select scans scan=trial['KcrsScans'], gaintype='KCROSS', solint='inf', refant=refant, smodel=[1,0,1,0], # Now we assume source IS polarized (doesn't matter what value) gaintable=[msdata+'.Bpass',msdata+'.G1'], interp=['nearest','linear']) # Apply so far (just for check) applycal(vis=msdata, calwt=False, gaintable=[msdata+'.Bpass',msdata+'.G1',msdata+'.Kcrs'], interp=['nearest','linear','nearest']) ################################### # XY phase calibration # ################################### # XY-phase and QU # (output caltable called 'XY0amb' because it may contain ambiguities) os.system('rm -rf '+msdata+'.XY0amb') gaincal(vis=msdata, caltable=msdata+'.XY0amb', gaintype='XYf+QU', # Solves for XYphase and Q & U (both % and PA) solint='inf', combine='scan,obs', # Combines 3 separate observations preavg=120, # In seconds -- a minimum parang change to detect polarization (?) refant=refant, smodel=[1,0,1,0], gaintable=[msdata+'.Bpass',msdata+'.G1',msdata+'.Kcrs'], interp=['nearest','linear','nearest']) # EXAMPLE TEXT: # Output from gaincal: # Spw = 0 (ich=32/64): # X-Y phase = 97.2269 deg. # Fractional Poln: Q = -0.0377898, U = -0.114666; P = 0.120733, X = -54.1202deg. # Net (over baselines) instrumental polarization: -0.00194052 # Use QU guess obtained from G1 above to fix ambiguities in XY0amb # (will use S for smodel below) os.system('rm -rf '+msdata+'.XY0') S=xyamb(xytab=msdata+'.XY0amb',qu=qu[0],xyout=msdata+'.XY0') # S is the correct model! Instead of, e.g., smodel=[1,0,1,0] # EXAMPLE TEXT: # Expected QU = (0.16099049532555496, 0.070475258915619607) # Spw = 0: Found QU = [-0.03778978 -0.11466616] # ...CONVERTING X-Y phase from 93.9271755416 to -86.0728244584 deg # Ambiguity resolved (spw mean): Q= 0.0377897769213 U= 0.114666156471 (rms= 0.0 0.0 ) P= 0.120732740711 X= 35.8798339194 # Returning the following Stokes vector: [1.0, 0.037789776921272278, 0.1146661564707756, 0.0] ########################################################################################## # Revise gain w/ good source pol estimate # # (cross-hands not used here, so cross-hand cal tables not needed) # # (NB: we are using S from above in smodel; parang=True so p-hand model has pol variation) # ########################################################################################## os.system('rm -rf '+msdata+'.G2') gaincal(vis=msdata, caltable=msdata+'.G2', # Now .G2 won't have polarization in gains solint='int', refant=refant, smodel=S, # Correct pol model gaintable=[msdata+'.Bpass'], interp=['nearest'], parang=True) # Critical for polarization calibration -- do in gaincal and applycal when applying polarization corrections -- .G2, .XY0 (XYphase), and .Df0 (D-terms) ###################################### # Apply so far, just to check # ###################################### applycal(vis=msdata, calwt=False, gaintable=[msdata+'.Bpass',msdata+'.G2',msdata+'.Kcrs',msdata+'.XY0'], # Took OUT .G1! Because now we have .G2 interp=['nearest','linear','nearest','nearest'], parang=True) ###################################### # D-term/leakage calibration # ###################################### os.system('rm -rf '+msdata+'.Df0') polcal(vis=msdata, caltable=msdata+'.Df0', solint='inf', combine='scan,obs', preavg=60, poltype='Dflls', smodel=S, refant=refant, # Solve for relative D-terms, because of discontinuity between spws (not an issue for calib) gaintable=[msdata+'.Bpass',msdata+'.G2',msdata+'.Kcrs',msdata+'.XY0'], interp=['nearest','linear','nearest','nearest']) ######################################## # Modify the D-term table in order for it # to be applied to the parallel-hands # as well as the cross-hands. ######################################## os.system('rm -rf '+msdata+'.Df0gen') Dgen(dtab=msdata+'.Df0',dout=msdata+'.Df0gen') ######################################### # Application of all calibration tables # ######################################### applycal(vis=msdata, calwt=False, gaintable=[msdata+'.Bpass',msdata+'.G2',msdata+'.Kcrs',msdata+'.XY0',msdata+'.Df0gen'], interp=['nearest','linear','nearest','nearest','linear'], parang=True) # Critical! ################################### # Split off fully calibrated data # ################################### os.system('rm -rf '+msdata_cal) os.system('rm -rf '+msdata_cal+'.flagversions') split(vis=msdata, datacolumn='corrected', outputvis=msdata_cal) ########### # Imaging # ########### os.system('rm -rf {1}.B{0}.IQUV.clean.*'.format(trial['band'],trial['source'])) clean(vis=msdata, imagename='{1}.B{0}.IQUV.clean'.format(trial['band'],trial['source']), weighting='briggs', robust=2, cell=['0.1arcsec'], imsize=trial['imsize'], mask='circle[[{0}pix,{0}pix], 2.0arcsec]'.format(trial['imsize'][0]/2.), stokes='IQUV', psfmode = 'clarkstokes', # Finds CLEAN components in each Stokes separately niter=300, interactive=False) ################################################ # Make Stokes I image (only) for selfcal model # ################################################ os.system('rm -rf {1}.B{0}.I.clean.*'.format(trial['band'],trial['source'])) clean(vis=msdata_cal, imagename='{1}.B{0}.I.clean'.format(trial['band'],trial['source']), weighting='briggs', robust=2, cell=['0.1arcsec'], imsize=trial['imsize'], mask='circle[[{0}pix,{0}pix], 2.0arcsec]'.format(trial['imsize'][0]/2.), stokes='I', psfmode = 'clarkstokes', niter=300, interactive=False) ########### # Selfcal # ########### # Not really necessary for any bands...but good practice. # Gaincal will use the MODEL column automatically # Do phase-only gaincal os.system('rm -rf '+msdata_cal+'.selfcal1') gaincal(vis=msdata_cal, caltable=msdata_cal+'.selfcal1', # First phase-only selfcal solution solint='int', # Bright source --> shortest solution interval refant=refant, calmode='p', minblperant=4, minsnr=4) # Plot solutions if interact: plotcal(msdata_cal+'.selfcal1', 'time', 'phase') # Apply phase-only selfcal solution applycal(vis=msdata_cal, calwt=False, gaintable=[msdata_cal+'.selfcal1'], interp=['linear']) ############ # Re-image # ############ os.system('rm -rf {1}.B{0}.IQUV.clean.selfcal1.*'.format(trial['band'],trial['source'])) clean(vis=msdata_cal, imagename='{1}.B{0}.IQUV.clean.selfcal1'.format(trial['band'],trial['source']), weighting='briggs', robust=2, cell=['0.1arcsec'], imsize=trial['imsize'], mask='circle[[{0}pix,{0}pix], 2.0arcsec]'.format(trial['imsize'][0]/2.), stokes='IQUV', psfmode = 'clarkstokes', niter=1000, interactive=False) ################## # Image analysis # ################## for stokes in ['I','Q','U','V'] : imfit(imagename='{1}.B{0}.IQUV.clean.selfcal1.image'.format(trial['band'],trial['source']), box='100,100,200,200', stokes=stokes, logfile='{2}.B{0}.IQUV.clean.image.stokes{1}.imfit'.format(trial['band'],stokes,trial['source'])) ###################################### # Make FITS file for reference image # ###################################### exportfits(imagename='{1}.B{0}.IQUV.clean.selfcal1.image'.format(trial['band'],trial['source']), fitsimage='{1}.B{0}.IQUV.clean.selfcal1.image.fits'.format(trial['band'],trial['source']) ) # Return the correct source model return S def offAxisRedux_onAxisCal(trial, smodel) : if not os.path.exists('./image_offaxis') : os.system('mkdir image_offaxis') # Check directory if os.path.abspath('.').split('/')[-1] != trial['dir'] : raise Exception("You're not in the correct directory!") # On-axis data msdata='onaxis.ms' msdata_cal='onaxis.cal.ms' # Get refant refant = trial['refant'] # Get datasets datasets = trial['datasets'] # Full, concatenated dataset msdata_off='concat.ms.split' # Gain table with on-axis corrections gaintable=[msdata+'.Bpass', msdata+'.Kcrs', msdata+'.G2', msdata+'.XY0', msdata+'.Df0gen', msdata_cal+'.selfcal1'] interp = ['nearest','nearest','linear','nearest','linear','linear'] # When starting over, use backup os.system('rm -rf concat.ms.split') os.system('cp -r concat.ms.split.backup concat.ms.split') ########################### # Prior to application of on-axis caltables to off-axis fields, # we need to calibrate the temporal phase variation, # which cannot be calibrated by on-axis calibration # # By using gaincal on all the data (msdata_off file), gaincal solves for # the phase of each time step, which includes all of the offaxis data # # In short, this is phase-only selfcal of each off-axis point. ###################################### # Calibrate residual phase (selfcal) # ###################################### os.system('rm -rf '+msdata_off+'.G3ph') gaincal(vis=msdata_off, caltable=msdata_off+'.G3ph', solint='int', refant=refant, calmode='p', smodel=smodel, gaintable=gaintable) # Add off-axis phase variation to gain table gaintable.append(msdata_off+'.G3ph') interp.append('linear') ##################################### # Application of calibration tables # ##################################### # Apply all corrections to all data applycal(vis=msdata_off, gaintable=gaintable, interp=interp, calwt=False, flagbackup=True, parang=True,) ########### # Imaging # ########### # Numbers for scan selection, imaging, etc. based on 11x11 size if trial['size'] == '11x11' : start = 100 stop = 223 delta = stop-start for i in range(start,stop): j=i-100 scan='{0},{1},{2}'.format(i, i+delta, i+2*delta) # Snagging same scan from each of the 3 runs # Image off-axis data with only on-axis leakages applied image = './image_offaxis/'+trial['source']+'.B'+str(trial['band'])+'.IQUV.clean.F'+str(j)+'' print ('CLEANing Stokes IQUV image {0}.'.format(j)) # Create off-axis images os.system('rm -rf {0}.image {0}.model {0}.psf {0}.residual {0}.flux'.format(image)) clean(vis=msdata_off, imagename=image, scan=scan, weighting='briggs', robust=2, cell=['0.1arcsec'], imsize=trial['imsize'], mask='circle[[{0}pix,{0}pix],2.0arcsec]'.format(trial['imsize'][0]/2.), stokes='IQUV', psfmode = 'clarkstokes', niter=1000, interactive=False) def offAxisImageAnalysis(trial) : ################## # Image analysis # ################## # Here we save the *.imfit files from 'DA', 'DV', and 'all' -- but we overwrite the offset images each time (to save space) # Figure out numbers for scan selection, imaging, etc. based on 11x11 vs. 9x9 sizes if trial['size'] == '11x11' : start = 100 stop = 223 delta = stop-start elif trial['size'] == '9x9' : start = 100 stop = 183 delta = stop-start # Full, concatenated dataset msdata_off='concat.ms.split' # Gaussian fit for IQUV images if not os.path.exists('./image_offaxis/analysis') : os.system('mkdir ./image_offaxis/analysis') for i in range(0,delta): print ('Fitting Stokes IQUV in image {0}.'.format(i)) for stokes in ['I','Q','U','V'] : image = './image_offaxis/'+trial['source']+'.B'+str(trial['band'])+'.IQUV.clean.F'+str(i)+'.image' imfit_file='./image_offaxis/analysis/'+trial['source']+'.B'+str(trial['band'])+'.IQUV.clean.F'+str(i)+'.image.stokes{0}.imfit'.format(stokes) os.system('rm -rf {0}'.format(imfit_file)) imfit(imagename=image, box='100,100,200,200', stokes=stokes, logfile=imfit_file) # Extract peak intensity of Stokes I,Q,U,V for stokes in ['I','Q','U','V'] : peak_file='./image_offaxis/analysis/{0}peak.B{1}.txt'.format(stokes,trial['band']) os.system('rm -rf {0}'.format(peak_file)) for i in range (1,delta-1): # Skips on-axis first/last scans; takes central value from the middle scan fr = open('./image_offaxis/analysis/'+trial['source']+'.B'+str(trial['band'])+'.IQUV.clean.F'+str(i)+'.image.stokes{0}.imfit'.format(stokes)) lines = fr.readlines() fr.close() i=str(i) success=False for line in lines: if line.find("Peak") >= 0: fw = open(peak_file, "a") fw.write(i) fw.write('\t') fw.write(line[:-1]) fw.write('\n') fw.close() success=True if not success : fw = open(peak_file, "a") fw.write(i) fw.write('\t') fw.write(' --- Peak: 0 +/- 0 mJy/beam (FAILED FIT)') fw.write('\n') fw.close() ############################################### # Copy *peak*txt files into current directory # ############################################### os.system('rm -rf ./{0}'.format(peak_file.strip('./image_offaxis/analysis/'))) os.system('cp -r {0} .'.format(peak_file))