############## script tested for CASA 6.2.1.7 - pipeline ############# ### For running a selection of a few steps use: #my_steps=[0,1,2,3,4,5,6,7,8,9,10,11,12,13] ### For running a row of many steps use: # my_steps=range(0,8) # from 0 to 7 ################# REQUIRES ANALYSIS UTILITIES #################### ### to obtain analysis scripts from NRAO site https://casaguides.nrao.edu/index.php/Analysis_Utilities ### open terminal and use following 2 bash commands: # wget ftp://ftp.cv.nrao.edu/pub/casaguides/analysis_scripts.tar # tar xvf analysis_scripts.tar ################################################################## ### uncomment and set pathToAU, if you haven't set up an init.py with this #pathToAU='YOUR-PATH-TO-AU' # #import sys #sys.path.append(pathToAU) #import analysisUtils as aU my_vis = 'uid___A002_X8a5fcf_X125f.ms' if 0 in my_steps : ## Import the asdm in .ms os.system('rm -rf '+my_vis) # erase existing '*.ms' os.system('rm -rf '+my_vis+'.flag*') # erase existing '*.falgversions' importasdm(asdm='../raw/uid___A002_X8a5fcf_X125f.asdm.sdm', # location of the raw data vis=my_vis, # name of the MeasurementSet to be created from the raw data lazy=True, # use symbolic link between raw data and MeasurementSet asis='*', # copy all metadata available (antennae position, pointing information bdfflags=True) # keep all existing flags of the raw data ## create summary of this EB os.system('rm -rf '+my_vis+'.listobs') # erase existing LISTOBS listobs(vis=my_vis, # name of the MS to be summarized listfile=my_vis+'.listobs') # Filename containing the summary ## Visualize the array configuration at the time of the observations os.system('rm -rf '+my_vis+'.plotants.png') # erase existing plot plotants(vis=my_vis, # name of the MS figfile=my_vis+'.plotants.png') # Filename of the output ##################################### #### Summary of listobs/plotants #### ##################################### # # ALMA Band-3 observations at ~107GHz with 4 SpWs with 1920 channel each # # # CALIBRATE FLUX/AMPLI: J0423-013 (ID:1) SpW: 17 - 19 - 21 - 23 # CALIBRATE BANDPASS : J0423-0120 (ID:0) SpW: 17 - 19 - 21 - 23 # CALIBRATE PHASE : J0423-0120 (ID:0) SpW: 17 - 19 - 21 - 23 # TARGET : NGC1614 (ID:3) SpW: 17 - 19 - 21 - 23 # # Tsys : Scan 3 - 5 - 8 SpW: 9 - 11 - 13 - 15 # (J0423-0120; J0423-013; NGC1614) # # # WVR : all scans SpW: 0 # # # ARRAY CONFIGURATION : 35 Antennae; Longest Baseline ~1150m -> Synthesize Beam ~0.5arcsec # Shortest Baseline ~30m -> LAS ~19.17arcsec # # ##################################### ##################################### ##################################### if 1 in my_steps : ##################################### ### A Priori Calibration: Tsys ### ##################################### ## Generation of the Tsys caltable using the task GENCAL # It is good practice to erase existing 'caltable' that # the forthcoming task will create. This avoids problems # when re-running a given calibration step os.system('rm -rf '+my_vis+'.tsys') # erase existing TSYS CALTABLE # generate of Tsys caltable (naming convention name_ms.caltable_name) gencal(vis=my_vis, # Name of the MS caltable=my_vis+'.tsys', # Name of the CALTABLE output by GENCAL caltype='tsys') # Tell GENCAL what we want to do # Tsys measurements are known to be corrupted at the edges # of each SpWs. We must flag these channels in the Tsys CALTABLE # using the task FLAGDATA # # CASA convention: # - '0:20,1:110' --> SpW 0 channel 20 and SpW1 channel 110 # - '0;2:20;110' --> SpW 0 and 2 and channel 20 and 110 # - '0~2:20;110' --> SpW 0 to 2 channel 20 and 110 # - '0~2:20;110,3:50 --> SpW 0 to 2 channel 20 and 110 and SpW 3 channel 50 # - '0~2:20~110,3:50 --> SpW 0 to 2 channel 20 to 110 and SpW 3 channel 50 flagdata(vis=my_vis+'.tsys', # Name of the MS or CALTABLE to be flagged mode='manual', # spw='9;11;13;15:0~3;124~127', # SpW and Channel to be flagged flagbackup = False) # Do we want to keep a backup of the flag before the run # We inspect these Tsys CALTABLE using the task PLOTBANDPASS, looking # for not well behaved antennae/scans. In case we find some, we will need # to flag the corresponding visibilities using FLAGDATA # plot for each antenna and each SPW, Tsys vs. Freq, overlaying all tsys scan together # and add the atmospheric transmission os.system('rm -rf '+my_vis+'.tsys.plots') # delete old folder for plots os.system('mkdir '+my_vis+'.tsys.plots') # make new folder for plots aU.plotbandpass(caltable=my_vis+'.tsys', # name of the CALTABLE to use xaxis='freq', yaxis='tsys', antenna='', # Which antennae do you want to inspect? spw='', # Which SpW do you want to inspect? showatm=True, # Show the atmospheric transmission? overlay='time', # figfile=my_vis+'.tsys.plots/'+my_vis, # convention: put it in a folder name_ms.caltable_name/name_ms interactive=False) # Open all the files generated by PLOTBANDPASS # >> eog my_files* & # # Inspect these plots and look for peculiar antennae/spw/scan # # Did you find any? if yes, those have to be flagged with FLAGDATA # These flags can be added in step 4 if 2 in my_steps : ##################################### ### A Priori Calibration: WVR ### ##################################### ## Generation of the WVR caltable using the task WVRGCAL # erase existing WVR CALTABLE os.system('rm -rf '+my_vis+'.wvr') # generate WVR CALTABLE wvrgcal(vis=my_vis, # Name of the MS to be process caltable=my_vis+'.wvr', # Name of the CALTABLE output by WRVGCAL spw=[17,19,21,23], # SpW to which solutions should be applied afterward smooth='6.05s', # Smooth the calibration on the given timescale tie=['J0423-0120,ngc_1614'], # Prioritise tieing the phase of these sources as well as possible (typically phase cal and target) statsource='ngc_1614') # Compute the statistics (Phase RMS, Disc) on this source only # Check the output of WVRGCAL in the CASA logger # Check the PWV value # PWV = xx mm # # As written in Help(WVRGCAL), Disc > few hundreds # indicates that WVR corrections should not be applied # # Did you find such cases? # # # Antenna with peculiar values should be flagged in WVRGCAL (wvrflag=xx). # # Do you find such antenna ? # # if 3 in my_steps : ############################################# ### A Priori Calibration: Application ### ############################################# ## Application of the WVR and Tsys CALTABLE using ## the task APPLYCAL # Application of these corrections on the field x applycal(vis=my_vis, # Name of the MS to process field='0~1', # field ID to correct spw='17,19,21,23', # SpW to be corrected (i.e., science SpWs) gaintable= [my_vis+'.tsys',my_vis+'.wvr'], # name of CALTABLEs to use gainfield= ['0~1',''], # ['Apply the Tsys from field x', 'Apply the WVR from field x'] interp ='linear,linear', # 'Time-interpolation method,Frequency-interpolation-method' calwt= True, # Applycal will scale weights of each visibility by 1/sqrt(Tsys(i) * Tsys(j)) flagbackup = False, spwmap=[[0,1,2,3,4,5,6,7,8,9,9,11,11,13,13,15,15,9,9,11,11,13,13,15,15],[]]) #[[apply Tsys SpW-x solutions to SpW-0, ... , Apply Tsys SpW-x solutions to SpW-i, ... , Apply Tsys SpW-x solutions to SpW-N],[same for WVR]] # your spw setup is # spws = [0,1,2,3,4,5,6,7,8] # spwmap applies the solution from the spw specified at a given list index in spwmap # to the spw at same list index in the spws list # spwmap = [0,0,2,2,4,4,6,7,8] ## Application of these corrections on the field x applycal(vis=my_vis, # Name of the MS to process field='2~3', # field ID to correct spw='17,19,21,23', # SpW to be corrected (i.e., science SpWs) gaintable= [my_vis+'.tsys',my_vis+'.wvr'], # name of CALTABLE to use gainfield= ['2',''], # [Apply the Tsys from field 'x', Apply the WVr from field 'x'] interp ='linear,linear', # 'Time-interpolation method,Frequency-interpolation-method' calwt= True, # Applycal will scale weights of each visibility by 1/sqrt(Tsys(i) * Tsys(j)) flagbackup = False, spwmap=[[0,1,2,3,4,5,6,7,8,9,9,11,11,13,13,15,15,9,9,11,11,13,13,15,15],[]]) #[[apply Tsys SpW-x solutions to SpW-0, ... , Apply Tsys SpW-x solutions to SpW-i, ... , Apply Tsys SpW-x solutions to SpW-N ],[same for WVR]] if 4 in my_steps : ############################################# ### A Priori Calibration: Flagging ### ############################################# ## We flag things which not anymore needed ## using the task FLAGDATA ### optional but skip for now: probably see suspicious antennas in autocorr ### # # # plotms(vis=my_vis, # # xaxis='frequency', # # yaxis='amp', # # field='0~2', # # spw='9,11,13,15', # Tsys channels # uvrange='', # # antenna='*&&&', # auto-correlations # avgtime='360', # scan length # avgfield=True, # average over all fields # iteraxis='baseline', # iterate over baseline # coloraxis='corr', # colour of data points = correlation # showgui=True, # # plotfile='field-0-1-2-amp-freq.png', # # overwrite=True) # # # ############################################################################### # We flag auto-correlation visibilities (i.e., antenna-i correlated with antenna-i) flagdata(vis = my_vis, # Name of the MS mode = 'manual', # Mode to be used spw = '1~24', # everything but the WVR SpW autocorr = True, flagbackup = False) # We flag un-necessary scans (POINTING, SIDEBAND_RATIO, ATMOSPHERE) flagdata(vis = my_vis, # Name of the MS mode = 'manual', intent = '*POINTING*,*SIDEBAND_RATIO*,*ATMOSPHERE*', flagbackup = False) # We flag visibilities coming from antennae which are shadowed by an other one. flagdata(vis = my_vis, mode = 'shadow', flagbackup = False) ## Additional flagging up to the students: ## Found while looking at Tsys flagdata(vis = my_vis, mode = 'manual', spw = '1~24', antenna='DA50', correlation='YY', flagbackup = False) if 5 in my_steps : ################################################### ### SPLIT I: slim down data for calibration ### ################################################### ## We now extract the Tsys and WVR corrected data of the science channels and drop the rest. ## In this way, the make the data easier to process for the main calibration. os.system('rm -rf '+my_vis+'.split') # erase existing '*.ms.split' mstransform(vis = my_vis, # Name of the input MS outputvis = my_vis+'.split', # Name of the output MS datacolumn = 'corrected', # only data with the the Tsys and WVR tables applied spw = '17,19,21,23', # science SpW typically keepflags = True) # keep flagged rows in the output ## We generate a listobs of the splitted data set which will be more tidy and easier to skim than the inital listobs. os.system('rm -rf '+my_vis+'.split.listobs') # erase existing LISTOBS listobs(vis = my_vis+'.split', listfile = my_vis+'.split.listobs') ## Scan and field IDs are the same, SpW Ids are renamed (compare frequencies)! # ## SpW 17 -> 0 ## SpW 19 -> 1 ## SpW 21 -> 2 ## SpW 23 -> 3 # ############################################################################## if 6 in my_steps : ###################################################### ### Calibration: Inspection, Editing, Flagging ### ###################################################### ## We have a look into the data with the plotms inspection tool. ## Depending on the selection of parameters we are plotting against each other, ## we can find problematic antennas, baselines, scans, correlatios, etc that ## we will most likely flag after inspection. ## Additional flagging up to the Students ## e.g. EDGE CHANNELS necessary? ## Found while looking at Amp vs UVDist (AFTER calibration!) #flagdata(vis = my_vis+'.split', # mode = 'xx', # spw='xx', # antenna='xx', # flagbackup = False) # as long as nothing is happening: pass #comment out as soon as you have something to flag above if 7 in my_steps : ########################################################## ### Calibration: set model for flux calibrator ### ########################################################## ################# QSO #################### ## Put a model to the AMPLITUDE/FLUX calibrator ## First, go to https://almascience.eso.org/sc/ ## Enter source name, frequency range of the observation (+ (10-20) GHz or larger or all bands), date of observation +/- one month. ## (Get values from the listobs-file) ## Derive the spectral index (S ~ nu^spix) and extrapolate the flux to the center of your frequency setup ## ############ Planet/Satellite ########### ## Just set a planet catalogue, e.g. Butler-JPL-Horizons 2012' as a standard and the field ID ## ########################################################################################################## ## ## From the ALMA calibrator catalogue we found that ## At nu_1=103.5GHz --> fnu_1=1.61 Jy ## At nu_2=91.5GHz --> fnu_2=1.67 Jy ## ## and fnu = A * nu^(spix) ## ## Let's solve for A and spix at 106.6GHz (central frequency of our observation; see listobs) ## ## ## fnu_1/fnu_2 = (nu_1/nu_2)^spix ## log10(fnu_1 / fnu_2) = spix * log10(nu_1 / nu_2) ## spix = log10(fnu_1 / fnu_2) / log10(nu_1 / nu_2) = log10( 1.61 / 1.67 ) / log10( 103.5 / 91.5 ) ## spix = - 0.29 ## ## so ## ## if nu_3 = 106.6GHz ## ## fnu_3 / fnu_2 = (nu_3/nu_2)^spix ## fnu_3 = fnu_2 * (nu_3/nu_2)^spix = 1.67 * (106.6 / 91.5)^(-0.29) = 1.60 Jy ## ## ############################################################################################################ setjy(vis = my_vis+'.split', # Name of the input MS standard = 'manual', # we do it field = 'J0423-013', # Flux calibrator fluxdensity = [1.60, 0, 0, 0], # Specified flux density at 'reffreq' in Jy [I,Q,U,V] - we don't have polarization data spix = -0.29, # spectral index reffreq = '106.6GHz') # reference frequency (usually middle of your spw setup) ## Check out the 'model' column in plotms (amp vs. uvdist and vs. frequency) if you are curious. plotms(vis = my_vis+'.split', xaxis = 'uvdist', yaxis = 'amp', ydatacolumn = 'model', field = 'J0423-013', avgchannel = '9999', coloraxis = 'spw', plotfile = my_vis+'.split.setjy.fluxcal.png') if 8 in my_steps : ##################################### ### Calibration: Bandpass ### ##################################### ## Check out the '*plotants.png' again. Which antenna can be used as reference antenna? ## It has to be at the center (i.e. covers large range of baselines), stable (amp/phase vs. time/freq), and with little or no flagging. ## First, get rid of the phase variations to improve the S/N. ## Reminder: Phase variations are combarably fast - choose the time interval on which you want to correct your data accordingly. os.system('rm -rf '+my_vis+'.split.ap_pre_bandpass') # delete old table in case it exists already gaincal(vis = my_vis+'.split', # Name of the input MS caltable = my_vis+'.split.ap_pre_bandpass', # Name of the output caltable field = '0', # J0423-0120 # bandpass calibrator spw = '', #scan = '4,7,10', solint = 'int', # solution interval (time interval) set to integration time refant = 'DA54', # at the center (i.e. covers large range of baseline), stable, with little or no flagging calmode = 'p') # calibrate phase ## Check the results in the plots below by iterating through the plotter: ## Look for discontinuities in phase vs. time, rapid wrapping of phase, ## large scatter in any solution... os.system('rm -rf '+my_vis+'.split.ap_pre_bandpass.plots')# delete old folder for plots os.system('mkdir '+my_vis+'.split.ap_pre_bandpass.plots') # make new folder for plots ## little tweak to get all antenna names tb.open(my_vis+'.split/ANTENNA') antList = tb.getcol('NAME') tb.close() ## iterate plots over antenna names for i in antList: plotms(vis=my_vis+'.split.ap_pre_bandpass', # name of the CALTABLE to use xaxis='time', # yaxis='phase', # field='0', # calibrator antenna=i, # needed to plot per antenna plotrange=[0, 0, -180, 180], # plotfile=my_vis+'.split.ap_pre_bandpass.plots/'+my_vis+'.split.ap_pre_bandpass.field0.'+i+'.png', gridrows=4, gridcols=1, customsymbol= True, #symbolshape='pixel', #symbolsize=2, #showlegend=True, iteraxis='spw', coloraxis='corr', exprange='all', # make plot for every iteration overwrite=True) ## Now, the BANDPASS correction! ## What were our assumptions on the time dependency? os.system('rm -rf '+my_vis+'.split.bandpass') # delete old table in case it exists already bandpass(vis = my_vis+'.split', # Name of the input MS caltable = my_vis+'.split.bandpass', # Name of the output caltable field = '0', # J0423-0120 # bandpass calibrator #scan = '4,7,10', solint = 'inf', # solution interval (time interval) set to infinite = entire scan combine = 'scan', # extend 'infinite' beyond scan limits ---> one solution per observation and channel refant = 'DA54', # reference antenna: at the center (i.e. covers large range of baseline), stable, with little or no flagging solnorm = True, # normalize solutions - easier to study/compare then - check description of this parameter's behaviour ('help') if bandpass is not the first calibration step! bandtype = 'B', # calibrate bandpass #degamp = 3, #degphase = 3, gaintable = my_vis+'.split.ap_pre_bandpass') # apply phase corrections on the fly before calculating bandpass corrections ## Check the results in the plots (folder) below: ## Look for discontinuities in phase vs. time, rapid wrapping of phase in bandpass solution, ## large roll-off in amplitude response near band edges in bandpass solution, large scatter in any solution... os.system('rm -rf '+my_vis+'.split.bandpass.plots') # delete old folder for plots os.system('mkdir '+my_vis+'.split.bandpass.plots') # make new folder for plots aU.plotbandpass(caltable=my_vis+'.split.bandpass', # name of the CALTABLE to use xaxis='freq', # what dependency did we correct for? yaxis='both', # amplitude and phase antenna='', # Which antenna to plot ? spw='', # showatm=True, # show atmosphere transmission curve field='0', # bandpass calibrator figfile=my_vis+'.split.bandpass.plots/'+my_vis+'.split.bandpass.field0', # convention: put it in a folder name_ms.caltable_name/name_ms interactive=False) if 9 in my_steps : ########################################### ### Calibration: Phase, Amplitude ### ########################################### ######### PHASE ######### ## First, find corrections for the phase variations on shortest timescales for all calibrators ## (The flux calibrator might require special treatment in case of planets and low signal/noise - beyond scope here). os.system('rm -rf '+my_vis+'.split.phase_int') # delete old table in case it exists already gaincal(vis = my_vis+'.split', # Name of the input MS caltable = my_vis+'.split.phase_int', # Name of the output caltable field = '0~1', # J0423-0120,J0423-013 # all calibrators solint = 'int', # solving for fast variations refant = 'DA54', # reference antenna gaintype = 'G', # gains for each polarization and spectral window calmode = 'p', # what quantitiy to calibrate gaintable = my_vis+'.split.bandpass') # apply bandpass corrections on the fly before calculating phase corrections for phase calibrator ## Check the results in the plots below by iterating through the plotter: ## Look for discontinuities in phase vs. time, rapid wrapping of phase, ## large scatter in any solution... os.system('rm -rf '+my_vis+'.split.phase_int.plots') # delete old folder for plots os.system('mkdir '+my_vis+'.split.phase_int.plots') # make new folder for plots ## little tweak to get all antenna names tb.open(my_vis+'.split/ANTENNA') antList = tb.getcol('NAME') tb.close() ## iterate plots over antenna names for i in antList: plotms(vis=my_vis+'.split.phase_int', # name of the CALTABLE to use xaxis='time', # yaxis='phase', # field='0,1', # calibrator antenna=i, # needed to plot per antenna plotrange=[0, 0, -180, 180], # plotfile=my_vis+'.split.phase_int.plots/'+my_vis+'.split.phase_int.field01.'+i+'.png', gridrows=4, gridcols=1, customsymbol= True, #symbolshape='circle', #symbolsize=1, title=i+' Phase Gain (int) vs. Time for field 0,1', #showlegend=True, iteraxis='spw', coloraxis='corr', exprange='all', # make plot for every iteration overwrite=True) ######### AMPLITUDE ######### ## Second, find corrections for the amplitude variations (slower compared to phase) ## on scan-long timescales (for good signal/noise) for all calibrators. os.system('rm -rf '+my_vis+'.split.ampli_inf') # delete old table in case it exists already gaincal(vis = my_vis+'.split', # Name of the input MS caltable = my_vis+'.split.ampli_inf', # Name of the output caltable field = '0~1', # J0423-0120,J0423-013 # all calibrators solint = 'inf', # solving for slower, scan block long variations refant = 'DA54', # reference antenna gaintype = 'T', # one gain solution for all polarizations calmode = 'a', # what quantitiy to calibrate gaintable = [my_vis+'.split.bandpass', my_vis+'.split.phase_int']) # apply bandpass and phase corrections on the fly before calculating amp corrections for phase calibrator ## Usually, at least for ALMA data, there is not much going on there, which is why ## we skip it during the lecture. ## Check the results in the plots below by iterating through the plotter: ## Look for discontinuities and large scatter in any solution... os.system('rm -rf '+my_vis+'.split.ampli_inf.plots') # delete old folder for plots os.system('mkdir '+my_vis+'.split.ampli_inf.plots') # make new folder for plots ## little tweak to get all antenna names tb.open(my_vis+'.split/ANTENNA') antList = tb.getcol('NAME') tb.close() ## iterate plots over antenna names for i in antList: plotms(vis=my_vis+'.split.ampli_inf', # name of the CALTABLE to use xaxis='time', # yaxis='amp', # field='0,1', # calibrator antenna=i, # needed to plot per antenna plotfile=my_vis+'.split.ampli_inf.plots/'+my_vis+'.split.ampli_inf.field01.'+i+'.png', gridrows=4, gridcols=1, customsymbol= True, #symbolshape='pixel', #symbolsize=2, title=i+' Amp Gain (inf) vs. Time for field 0,1', #showlegend=True, iteraxis='spw', coloraxis='corr', exprange='all', # make plot for every iteration overwrite=True) ######### PHASE again! ######### ## Third, find corrections for the phase variations on scan-long timescales for all calibrators or at least the phase calbrator. ## We need this table for the target source. High time resolution (solint='int') is useless. ## Instead we want an average value for each phase calibrator scan before and after the target scan. os.system('rm -rf '+my_vis+'.split.phase_inf') # delete old table in case it exists already gaincal(vis = my_vis+'.split', # Name of the input MS caltable = my_vis+'.split.phase_inf', # Name of the output caltable field = '0~1', # J0423-0120,J0423-013 # all calibrators solint = 'inf', # solving for slower, scan block long variations refant = 'DA54', # reference antenna gaintype = 'G', # gains for each polarization and spectral window calmode = 'p', # what quantitiy to calibrate gaintable = my_vis+'.split.bandpass') # apply bandpass corrections on the fly before calculating phase corrections for phase calibrator ## Usually, at least for ALMA data, there is not much going on there, which is why ## we skip it during the lecture. ## Check the results in the plots below by iterating through the plotter: ## Look for discontinuities in phase vs. time, rapid wrapping of phase, ## large scatter in any solution... os.system('rm -rf '+my_vis+'.split.phase_inf.plots') # delete old folder for plots os.system('mkdir '+my_vis+'.split.phase_inf.plots') # make new folder for plots ## little tweak to get all antenna names tb.open(my_vis+'.split/ANTENNA') antList = tb.getcol('NAME') tb.close() ## iterate plots over antenna names for i in antList: plotms(vis=my_vis+'.split.phase_inf', # name of the CALTABLE to use xaxis='time', # yaxis='phase', # field='0,1', # calibrator antenna=i, # needed to plot per antenna plotrange=[0, 0, -180, 180], # plotfile=my_vis+'.split.phase_inf.plots/'+my_vis+'.split.phase_inf.field01.'+i+'.png', gridrows=4, gridcols=1, customsymbol= True, #symbolshape='pixel', #symbolsize=2, title=i+' Phase Gain (inf) vs. Time for field 0,1', #showlegend=True, iteraxis='spw', coloraxis='corr', exprange='all', # make plot for every iteration overwrite=True) if 10 in my_steps : ########################################## ### Calibration: FLUX TRANSFER ### ########################################## ## rescale the amplitudes to real fluxes based on the model for your flux calibrator. os.system('rm -rf '+my_vis+'.split.flux_inf') # delete old table in case it exists already fluxscale(vis = my_vis+'.split', # Name of the input MS caltable = my_vis+'.split.ampli_inf', # Name of the input caltable fluxtable = my_vis+'.split.flux_inf', # Name of the output caltable reference = '1') # J0423-013 # flux calibrator - field from which the fluxscale is transferred ## Double check the resulting fluxes for the other (non-flux) calibrators in the CASA logger ## with the observed fluxes from the catalogue https://almascience.eso.org/sc/. if 11 in my_steps : ############################################################ ### Calibration: Apply the gaintables to the data ### ############################################################ ## Apply the solutions for the flux calibrator to its data applycal(vis = my_vis+'.split', # Name of the input MS field = '1', #J0423-013 FLUX # Flux calibrator gaintable = [my_vis+'.split.bandpass', my_vis+'.split.phase_int', my_vis+'.split.flux_inf'], # apply bandpass, fast phase and flux correction gainfield = ['', '1', '1'], # use corrections from the calibrator in 'field', except from the bandpass table - wwhat to set here? interp = 'linear,linear', calwt = True, # Calibrate weights along with data for each gaintable flagbackup = False) ## Apply the solutions for the phase calibrator to its data (we can include the bandpass calibrator here since it is the same source). applycal(vis = my_vis+'.split', # Name of the input MS field = '0', # J0423-0120 BANDPASS and PHASE # phase (and bandpass) calibrator; bandpass would need an own applycal command if fields were not the same gaintable = [my_vis+'.split.bandpass', my_vis+'.split.phase_int', my_vis+'.split.flux_inf'], # apply bandpass, fast phase and flux correction gainfield = ['', '0', '0'], # use corrections from the calibrator in 'field', except from the bandpass table - wwhat to set here? interp = 'linear,linear', calwt = True, # Calibrate weights along with data for each gaintable flagbackup = False) ## Apply the solutions for the phase calibrator to the target data applycal(vis = my_vis+'.split', # Name of the input MS field = '3', # ngc_1614 TARGET # target source gaintable = [my_vis+'.split.bandpass', my_vis+'.split.phase_inf', my_vis+'.split.flux_inf'], # apply bandpass, slow phase and flux correction gainfield = ['', '0', '0'], # use corrections from the PHASE calibrator ! , except from the bandpass table - wwhat to set here? interp = 'linear,linear', calwt = True, # Calibrate weights along with data for each gaintable flagbackup = False) if 12 in my_steps : ###################################################### ### Calibration: Post-Calibration Inspection ### ###################################################### ## We have a look into the CORRECTED data with the plotms inspection tool again. #os.system('rm -rf '+my_vis+'.split.inspect.plots') # delete old folder for plots - uncomment on purpose! os.system('mkdir '+my_vis+'.split.inspect.plots') # make new folder for plots - might contain images from several script calls flds=['0','1'] allcals = ','.join(flds) corrrun='uncorr' #corrrun='1st_corr' #corrrun='2nd_corr' datacol='data' #datacol='corrected' ## Amplitudes/Phase vs. UVDistances (repeat for all fields) for fld in flds: plotms(vis=my_vis+'.split', # Name of the input MS xaxis='uvdist', yaxis='amp', # redo for phase (interactively) ydatacolumn=datacol, # use the column of the corrected data field=fld, # repeat for every field, but at least calibrators (interactively) avgchannel='1e8', # averaging over ... avgtime='1e8', avgscan=True, iteraxis='spw', # plot for each SpW (iteration) - redo for iteration over scan (set also avgscan=False!) coloraxis='corr', # colour code for data points showgui=True, # open user interface plotfile=my_vis+'.split.inspect.plots/field-'+fld+'-amp-uvdist_'+corrrun+'.png', # export the plot - iteration name will be inserted automatically - change field name manually! exprange='all', # make plot for every iteration overwrite=True) plotms(vis=my_vis+'.split', # Name of the input MS xaxis='uvdist', yaxis='phase', # redo for phase (interactively) ydatacolumn=datacol, # use the column of the corrected data field=fld, # repeat for every field, but at least calibrators (interactively) avgchannel='1e8', # averaging over ... avgtime='1e8', avgscan=True, iteraxis='spw', # plot for each SpW (iteration) - redo for iteration over scan (set also avgscan=False!) coloraxis='corr', # colour code for data points showgui=True, # open user interface plotfile=my_vis+'.split.inspect.plots/field-'+fld+'-ph-uvdist_'+corrrun+'.png', # export the plot - iteration name will be inserted automatically - change field name manually! exprange='all', # make plot for every iteration overwrite=True) # Most data should fall along a straight line (for unresolved or marginally resolved sources) # or should follow something similar to a sinc^2-function (for solar system objects). # Look for outliers to check. ######## Note to lecturer ###### ## ## amplitude after calibration with no flags: only DA41 spw0 is out of line ## ## phase after calibration with no flags: between +/-2 deg ## ################################ ## Amplitudes/Phase vs. Channels (repeat for all fields) plotms(vis=my_vis+'.split', xaxis='frequency', yaxis='amp', # redo for phase (interactively) field=fld, # repeat for every field, but at least calibrators (interactively) ydatacolumn=datacol, # use the column of the corrected data avgtime='1e8', avgscan=True, iteraxis='spw', # plot for each SpW (iteration) - redo for iteration over scan (set also avgscan=False!) coloraxis='corr', showgui=True, plotfile=my_vis+'.split.inspect.plots/field-'+fld+'-amp-channel_'+corrrun+'.png', # export the plot - iteration name will be inserted automatically - change field name manually! overwrite=True) ## DV07 bad in spw0 f3 plotms(vis=my_vis+'.split', xaxis='frequency', yaxis='phase', # redo for phase (interactively) field=fld, # repeat for every field, but at least calibrators (interactively) ydatacolumn=datacol, # use the column of the corrected data avgtime='1e8', avgscan=True, iteraxis='spw', # plot for each SpW (iteration) - redo for iteration over scan (set also avgscan=False!) coloraxis='corr', showgui=True, plotfile=my_vis+'.split.inspect.plots/field-'+fld+'-ph-channel_'+corrrun+'.png', # export the plot - iteration name will be inserted automatically - change field name manually! overwrite=True) # Most data should fall along a relatively straight line or look smooth. # Channels at SpW edge can be unusually high or low (bandpass roll-off). # Look for roll-offs, bad data (outliers) and (potentially) atmospheric absorption features to check. ######## Note to lecturer ###### ## ## all smooth but minor offsets between correlations ## ################################ ## Amplitudes/Phase vs. Time (repeat for all fields) #corrrun='1st_corr' plotms(vis=my_vis+'.split', field='', xaxis='time', yaxis='amp', # redo for phase (interactively) ydatacolumn=datacol, # use the column of the corrected data avgchannel='1e8', iteraxis='spw', coloraxis='field', showgui=True, plotfile=my_vis+'.split.inspect.plots/field-cal-amp-time_'+corrrun+'_inclTarget.png', exprange='all', # make plot for every iteration overwrite=True) plotms(vis=my_vis+'.split', field=allcals, xaxis='time', yaxis='phase', # redo for phase (interactively) ydatacolumn=datacol, # use the column of the corrected data avgchannel='1e8', iteraxis='spw', coloraxis='field', showgui=True, plotfile=my_vis+'.split.inspect.plots/field-cal-ph-time_'+corrrun+'.png', exprange='all', # make plot for every iteration overwrite=True) # Amplitude: Each field should cover an amplitude range constant over entire observation. # Outliers are bad integrations/scans. # A gradual change could be related to the elevation (to check), else investigate. # Phase: The data should look smooth enough to interpolate between the phase calibration observations. # Gradual changes in phase are ok, but jumps (e.g. >120 deg) between phase calibration observations can be problematic. # Abnormal scatter indicates bad data. # Look for suspicious behaviour. ######## Note to lecturer ###### ## ## amplitudes after calibration with no flags: ints with low amp still present, even when !DA45;!DA64;!DA41;!DV14 ## ## phase after calibration with no flags: between +/-2 deg ## ################################ if 13 in my_steps : ########################################### ### SPLIT II: your calibrated data ! ### ########################################### ## We now extract the corrected the data column for all sources and carve them in stone. ## SPLIT, prepare data for imaging os.system('rm -rf '+my_vis+'.split.cal') # delete old table in case it exists already mstransform(vis = my_vis+'.split', # Name of the input MS outputvis = my_vis+'.split.cal', # Name of the output MS datacolumn = 'corrected', # only thge calibrated and data keepflags = True) # keep flagged rows in the output ## Congrats to your calibrated and reduced data set! We can now move on to IMAGING os.system('rm -rf '+my_vis+'.split.cal.listobs') # erase existing LISTOBS listobs(vis = my_vis+'.split.cal', listfile = my_vis+'.split.cal.listobs')