# ALMA Solar Science-Data Reduction Script # # $Id: uid.ms.scriptForCalibration.py,v 1.1 2019/10/25 12:18:50 dpetry Exp $ # # Original template by Masumi Shimojo (EA ARC) # Modified by Miro Barta (EU ARC, Czech node) on 2018/07/21 # Modified by masumi.shimojo@nao.ac.jp (EA ARC) on 2018/07/11 for CASA 5.3.0 # Step 20 changed by barta@asu.cas.cz (EU ARC) on 2018/07/28 # Modified for CASA 5.4-70 by barta@asu.cas.cz (EU ARC) on 2019/04/21 # Modified for CASA 5.6.1-8 by barta@asu.cas.cz (EU ARC) on 2019/11/21 # Modified by masumi.shimojo@nao.ac.jp (EA ARC) on 2019/10/02 for Band 7 # Revised for CASA6+/Python3 by barta@asu.cas.cz on 2023/12/27 # Updated by Tant plots by barta@asu.cas.cz on 2025/12/12 # # Working in CASA version: 6.0.0 and higher # # #========================================================================= # # !!! Parts that are of special attention of the data analysts # are marked by _SADA_ string in the comments !!! # #========================================================================= # Initial actions #================================================== # Switch for each process thesteps = [] step_title = {0: 'Import of the ASDM', 1: 'Print observation summary info to files', 2: 'Fix of SYSCAL table times', 3: 'A priori flagging1', 4: 'Generation of the Tsys cal tables', 5: 'Generation of antenna-position correction table', 6: 'Application of Tsys, AntPos, & Tant to calibrators and Sun', 7: 'A priori flagging2', 8: 'Split out science SPWs', 9: 'Initial flagging', 10: 'Save original flags', 11: 'Putting a model for the flux calibrator', 12: 'Save flags - restart point for calibration', 13: 'Bandpass calibration', 14: 'Save flags before gain cal', 15: 'Gain calibration', 16: 'Save flags before applycal', 17: 'Application of the bandpass and gain cal tables', 18: 'Split out corrected column', 19: 'Save flags after applycal', 20: 'Re-calculation of the direction'} #================================================== # The Python variable 'mysteps' will control which steps # are executed when you start the script using # execfile('scriptForCalibration.py') # e.g. setting # mysteps = [2,3,4]# before starting the script will make the script execute # only steps 2, 3, and 4 # Setting mysteps = [] will make it execute all steps. try: print('List of steps to be executed ...', mysteps) thesteps = mysteps except: print ('global variable mysteps not set.') if (thesteps==[]): thesteps = range(0,len(step_title)) print ('Executing all steps: ', thesteps) #================================================== # Check CASA ver. import casalith cversion = casalith.version_string() print ("You are using " + cversion) if (cversion < '6.0.0'): print ("YOUR VERSION OF CASA IS TOO OLD FOR EXECUTING THIS SCRIPT.") print ("PLEASE UPDATE IT BEFORE PROCEEDING.") else: print ("Your CASA is appropriate for data processing using this script.") #================================================== # Parameters Setup #--- ASDM Directory Name __KEY_ASDM_NAME__ # SSG will replace this key with actual line(s) of code! #--- Measurement Sets names # Observed MS mso = asdm + '.ms' # MS after split-out of science SPWs mss = asdm + '.ms.split' # Calibrated MS msc = mss + '.cal' #--- Reference antenna # Reference antenna as auto-detected by es.getRefAntenna(msName) __KEY_REF_ANT__ # SSG will replace this key with actual line(s) of code! # _SADA_ # # It may happen that auto-detected reference antenna above has problems (e.g., # many or long down-time breaks, high TSys or XX/YY TSys disbalance, etc). # In such a case, choose the antena manually: Select one with little breaks, # close to the center, with good Tsys && Tant - use .flagcmd.png plot, # plotants() and Tsys, Tsys.ext, and Tant plots. # Commented example: # refAnt = 'DV23' #--- Fields/scans for calibrators and the science target. Autodetected using # msmd.fieldsforintent() and msmd.scansforintent() __KEY_OBS_FIELDS__ # SSG will replace this key with actual line(s) of code! #================================================== # Loading general-use modules. import re import os import csv #================================================== # Loading the solar-specific utilities. import glob # This is the first guess of their location (local, in the project dir) sunUtils = glob.glob('../sun_reduction_util.py') if len(sunUtils) == 0: # Not found locally. # Search in the 'script' subdir; let get prepended the specific MOUS id. sunUtils = glob.glob('../../script/*sun_reduction_util.py') execfile(sunUtils[0],globals()) #========================================================================= # Calibration step-by-step. #========================================================================= # Step 0: Import of the ASDM mystep = 0 if(mystep in thesteps): casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO') print ('Step ', mystep, step_title[mystep]) if not os.path.exists(mso): importasdm(asdm = asdm, vis = mso, lazy = True, asis='Antenna Station Receiver Source CalAtmosphere ' +'CalWVR CorrelatorMode SBSummary CalDevice') else: print("MS has been already imported. Delete the MS dir if re-import " \ +"is required.") #========================================================================= # Step 1: Create the files with observation summary mystep = 1 if(mystep in thesteps): casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO') print('Step ', mystep, step_title[mystep]) # Dump of the observing information: CASA logger + file. os.system('rm -rf '+ asdm + '_listobs.txt') listobs(mso, listfile = asdm + '_listobs.txt') __KEY_LISTOBS_PLOT__ # es.listobs3() in QA2 mode, empty otherwise. #========================================================================= # Step 2: Fix of SYSCAL table times mystep = 2 if(mystep in thesteps): casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO') print ('Step ', mystep, step_title[mystep]) from casarecipes.almahelpers import fixsyscaltimes fixsyscaltimes(vis = mso) #========================================================================= # Step 3: A priori flagging1 mystep = 3 if(mystep in thesteps): casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO') print('Step ', mystep, step_title[mystep]) # Data with the intent that is no more required flagdata(vis = mso, mode = 'manual', intent = '*POINTING*,*SIDEBAND_RATIO*', flagbackup = False) # Channel-integrated SPWs 5,7,9,11 flagdata(vis = mso, mode = 'manual', spw = '6,8,10,12', flagbackup = False) # Show (and possibly apply) flags for antenna downtimes. flagcmd(vis = mso, inpmode = 'table', useapplied = True, action = 'plot', plotfile = mso+'.flagcmd.png') ''' # This may lead to problems with lost of too much data - be carefull using it. flagcmd(vis = mso, inpmode = 'table', useapplied = True, action = 'apply') ''' # Save so-far done flags flagmanager(vis = mso, mode = 'save', versionname = 'priori1') #========================================================================= # Step 4: Generate the Tsys caltable mystep = 4 if(mystep in thesteps): casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO') print('Step ', mystep, step_title[mystep]) os.system('rm -rf ' + mso + '.tsys') os.system('rm -rf ' + mso + '.tsys_org') # Generate Tsys caltable gencal(vis = mso, caltable = mso + '.tsys', caltype = 'tsys') # Flagging edge channels in Tsys table flagdata(vis = mso + '.tsys', mode = 'manual', spw = '5:0~9;116~127,7:0~9;116~127,9:0~9;116~127,11:0~9;116~127', flagbackup = False) __KEY_TSYS_PLOT__ # SSG: Tsys plots in the QA2 regime, empty otherwise. #========================================================================= # Step 5: Generation of the antenna position cal table mystep = 5 if(mystep in thesteps): casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO') print('Step ', mystep, step_title[mystep]) # Remove an old caltable, if it exists. os.system('rm -rf '+mso+'.antpos') # Generate the new antena-position correction cal table. __KEY_ANT_POS_TABLE__ # SSG will replace this key with actual line(s) of code! #========================================================================= # Step 6: Application of Tsys and AntPos tables to the calibrators. # Generation of T_ant gaintables and application of Tsys, Tant, # and antPos corrections to the Sun. mystep = 6 if(mystep in thesteps): casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO') print ('Step ', mystep, step_title[mystep]) from casarecipes.almahelpers import tsysspwmap tsysmap = tsysspwmap(vis = mso, tsystable = mso+'.tsys', tsysChanTol = 1) #------- Apply Tsys and TantPos corrections to all calibrators # # Note: We are using gainfield='nearest' in order to handle all calibrators # homogeneously, either they have their own Tsys scans or not. __KEY_APPLY_PRECALIBS__ # SSG replaces this key with actual line(s) of code #------- Do T_ant calculations for the Sun (solar ObsMode specific). # # Explanatory notes: # sol_ampcal_2 from the sun_reduction_utils.py has a tripple effect: # (1) It can calculate the Tant caltable (i.e. doing sort of gaincal()) # (2) It applies the Tant caltable generated by itself and the Tsys & antPos # caltables created in Steps 4 & 5, respectivelly ( call to applycal()). # (3) It can generate CSV files with antenna temperatures Tant that are # usefull for diagnostic plots via sol_tsunplot(asdm). # # The generation of Tant caltable is quite time-demanding. Hence, if the # CSV files and the Tant caltable already exists, the sol_ampcal_2() can be # called with parameters # sol_ampcal_2(exisTbl=True, outCSV=False) # in order to call just the applycal(). # Call with gaincal(), CSV generation and applycal() sol_ampcal_2(mso, mso+'.tsys', mso+'.antpos',exisTbl=False, outCSV=True) # Call with applycal() only. Use only if you _do_ have the T_ant gaintables # from your previous run (i.e. you are repeating this step for some reason). # sol_ampcal_2(mso, mso+'.tsys', mso+'.antpos', exisTbl=True, outCSV=False) __KEY_CAL_TANT_PLOT__ # Tant and calWeight plots for QA2, empty otherwise. #========================================================================= # Step 7: A priori flagging2 mystep = 7 if(mystep in thesteps): casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO') print('Step ', mystep, step_title[mystep]) # Tsys calibration is over - we do not need ATM data flagdata(vis = mso, mode = 'manual', intent = '*ATMOSPHERE*', flagbackup = False) # Flag autocorrelations flagdata(vis = mso, mode = 'manual', autocorr = True, flagbackup = False) # WVR data - not to be used for solar observations flagdata(vis = mso, mode = 'manual', spw = '0~4', flagbackup = False) # Other now useles SPWs - the science SPWs are 5,7,9, and 11 flagdata(vis = mso, mode = 'manual', spw = '13~36', flagbackup = False) # Flagging edge channels (TDM) flagdata(vis = mso, spw='5:0~9;116~127,7:0~9;116~127,9:0~9;116~127,11:0~9;116~127', mode = 'manual', flagbackup = False) # Flagging shadowed data. # Sun elevation limits (40 - 70 deg) actually most likely prevents shadowing. # flagdata(vis = mso, mode = 'shadow', flagbackup = False) # _SADA_ #--- Flagging antennas with bad/missing Tsys or Tant: # To be modified by the analyst. The code here is just an example: ''' # CM01: Too low T_ant for BB=3,4 -> SPW=9,11; XX pol.; scan 9 flagdata(vis = mso, antenna = 'CM01', spw = '9,11', correlation='XX', scan='9', mode = 'manual', flagbackup = False) ''' # Save flags done so far. flagmanager(vis = mso, mode = 'save', versionname = 'priori2') #========================================================================= # Step 8: Split out science SPWs mystep = 8 if(mystep in thesteps): casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO') print('Step ', mystep, step_title[mystep]) os.system('rm -rf ' + mss) os.system('rm -rf '+ mss + '.flagversions') # Keep only data with usefull intent. # NB: The P_off solar sub-scans are dropped here by ON_SOURCE specification __KEY_SPLIT_INTENT__ # SSG replaces this key with actual line(s) of code. split(vis = mso, outputvis = mss, datacolumn = 'corrected', intent=sciScansIntent, spw = '5,7,9,11', keepflags = True) #========================================================================= # Step 9: Initial flagging mystep = 9 if(mystep in thesteps): casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO') print('Step ', mystep, step_title[mystep]) # Flagging zero values flagdata(vis = mss, mode = 'clip', clipzeros = True, flagbackup = False) # Flagging Earth-atmosferic lines (if there are some): __KEY_FLAG_ATM_LINES__ # SSG replaces this key with actual line(s) of code #========================================================================= # Step 10: Save original flags mystep = 10 if(mystep in thesteps): casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO') print('Step ', mystep, step_title[mystep]) if not os.path.exists(mss + '.flagversions/Original.flags'): flagmanager(vis = mss, mode = 'save', versionname = 'Original') #========================================================================= # Step 11: Putting a model for the flux calibrator(s) mystep = 11 if(mystep in thesteps): casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO') print ('Step ', mystep, step_title[mystep]) # Set the flux-density (and spec-slope) for the main flux calibrator. # It has been found via aU.getALMAFlux() __KEY_SET_JY__ # SSG will replace this key with actual line(s) of code! #========================================================================= # Step 12: Save flags - restart point before calibration mystep = 12 if(mystep in thesteps): casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO') print ('Step ', mystep, step_title[mystep]) # _SADA_ # Here the analyst can put flagging of bad antennas as found by inspection # of G-cal/B-cal plots created in subsequent steps 13-15. # If additional flagging is set here, the calibration procedure can be # restarted/iterated from this point this way # in CASA: # flagmanager(vis = mss, mode = 'restore', versionname = 'Original') # mysteps=range(12,16) # execfile('uid__.....ms.scriptForCalibration.py') ''' # This is just an example !!! # # Wrong already pre-bandpass caltable (and consequently all following steps) # for antenna CM02, spw=1 flagdata(vis = mss, antenna = 'CM02', spw = '1', mode = 'manual', flagbackup = False) ''' # Save the flags done so far flagmanager(vis = mss, mode = 'save', versionname = 'BeforeBandpassCalibration') #========================================================================= # Step 13: Bandpass calibration mystep = 13 if(mystep in thesteps): casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO') print('Step ', mystep, step_title[mystep]) os.system('rm -rf '+ mss +'.ap_pre_bandpass*') # Calibration of fast phase distortions over the scan of band-pass calibrator gaincal(vis = mss, caltable = mss + '.ap_pre_bandpass', field = bpassCalField, scan = bpassCalScan, solint = 'int', refant = refAnt, calmode = 'p') os.system('rm -rf '+mss+'.bandpass*') # Band-pass calibration. Apply correction of fast phase variations (done # above) on-the-fly using the gaintable parameter. bandpass(vis = mss, caltable = mss+'.bandpass', field = bpassCalField, scan = bpassCalScan, solint = 'inf', refant = refAnt, solnorm = True, bandtype = 'B', gaintable = mss+'.ap_pre_bandpass') __KEY_BP_CAL_PLOT__ # SSG: Make BP-cal diagnostic plots in the QA2 regime. #========================================================================= # Step 14: Save flags before gain cal mystep = 14 if(mystep in thesteps): casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO') print('Step ', mystep, step_title[mystep]) flagmanager(vis = mss, mode = 'save', versionname = 'BeforeGainCalibration') #========================================================================= # Step 15: Gain calibration mystep = 15 if(mystep in thesteps): casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO') print('Step ', mystep, step_title[mystep]) # barta@asu.cas.cz: This step is patterned according to CASA guide at # https://casaguides.nrao.edu/index.php/AntennaeBand7_Calibration # In this respect it differs from what the Eric Villard's script generator # has put here. # In particular, gaintype='G' (_not_ 'T') and calmode='ap' (not 'a') is used # for amplitude caltable on the scan-time scale. See also Masumi's analysis # of XX/YY asymmetric distributions and 'T' vs. 'G' gaintype: # https://www.evernote.com/l/AKegK4x_iKtKkL3pqYQK98q1WUPgK-5msEE. #--- Fast (time scale of integration/subscan) phase variations # Apply band-pass corrections on-the_fly os.system('rm -rf '+mss+'.phase_int*') gaincal(vis = mss, caltable = mss + '.phase_int', field = allCals, solint = '1.92s', # todo: Find it by the SSG refant = refAnt, gaintype = 'G', calmode = 'p', minsnr = 2.0, # todo: Find-out by the SSG gaintable = mss + '.bandpass') #--- Slow (on the scan timescale) phase variations. # Apply band-pass corrrections on-the_fly os.system('rm -rf '+mss+'.phase_inf*') gaincal(vis = mss, caltable = mss+'.phase_inf', field = allCals, solint = 'inf', preavg=1.92, # todo: Find it by the SSG refant = refAnt, gaintype = 'G', calmode = 'p', minsnr = 3.0, gaintable = mss+'.bandpass') #--- Slow (timescale of scan) amplitude variation. # Apply the bandpass and fast phase variations on-the-fly os.system('rm -rf '+ mss +'.ampli_inf*') gaincal(vis = mss, caltable = mss + '.ampli_inf', field = allCals, solint = 'inf', preavg=1.92, # todo: Find it by the SSG refant = refAnt, gaintype = 'G', calmode = 'ap', minsnr = 3.0, gaintable = [mss + '.bandpass', mss + '.phase_int']) #--- Set absolute scale to all calibrators os.system('rm -rf '+mss+'.flux_inf*') os.system('rm -rf '+mss+'.fluxscale*') # Log temporarilly to the 'uid__.. ms.split.fluxscale' file mylogfile = casalog.logfile() casalog.setlogfile(mss + '.fluxscale') # Bootstrap the flux from Flux/Amp calibrator to other (P-cal, B-cal if it is # not the same). fluxscaleDict = fluxscale(vis = mss, caltable = mss + '.ampli_inf', fluxtable = mss + '.flux_inf', reference = ampCalField) # Reset logfile back to the default. casalog.setlogfile(mylogfile) __KEY_G_CAL_PLOT__ # SSG: Make gain-cal diagnostic plots in the QA2 regime. #========================================================================= # Step 16: Save flags before applycal mystep = 16 if(mystep in thesteps): casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO') print('Step ', mystep, step_title[mystep]) flagmanager(vis = mss, mode = 'save', versionname = 'BeforeApplycal') #========================================================================= # Step 17: Application of the bandpass and gain cal tables mystep = 17 if(mystep in thesteps): casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO') print('Step ', mystep, step_title[mystep]) #--- Apply the cal tables to F-cal (if we have one separate), B-cal and P-cal for i in allCalsList: applycal(vis = mss, field = str(i), gaintable = [mss+'.bandpass',mss+'.phase_int',mss+'.flux_inf'], gainfield = [bpassCalField, i, i], interp=['nearest, linear','nearest','nearest'], calwt = True, flagbackup = False) #--- Apply calibrations to the science target # Use the phase calibrator as its observations interleave that of the Sun applycal(vis = mss, field = targetFields, # See definition in the script header. gaintable = [mss+'.bandpass', mss+'.phase_inf', mss+'.flux_inf'], gainfield = [bpassCalField, phaseCalField, phaseCalField], interp=['nearest, linear','linear','linear'], calwt = True, flagbackup = False) #========================================================================= # Step 18: Split out corrected column mystep = 18 if(mystep in thesteps): casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO') print('Step ', mystep, step_title[mystep]) os.system('rm -rf '+msc+'*') os.system('rm -rf '+msc+'*.flagversions') split(vis = mss, outputvis = msc, datacolumn = 'corrected', keepflags = True) #========================================================================= # Step 19: Save flags after applycal mystep = 19 if(mystep in thesteps): casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO') print('Step ', mystep, step_title[mystep]) flagmanager(vis = msc, mode = 'save', versionname = 'AfterApplycal') #========================================================================= # Step 20: Setting sky-coordinates (RA/Dec) to the target fields. mystep = 20 if(mystep in thesteps): casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO') print('Step ', mystep, step_title[mystep]) #--- Use fixplanets() for attributing the RA/Dec coordinates to the solar # fields at the reference time. # Sky positions (RA/Dec) for individual solar fields. They have been found # via msmd.phasecenter() at reference time. __KEY_SUN_FIELD_DIRS__ # SSG will replace this key with actual lines of code! # Calculate the phase directions for all Sun fields. for theField in SunFields: # Set the phasecenters found for the field by msmd.phasecenter() to the # FIELD table via fixplanets() - overwrite the original OFFSETs fixplanets(msc,field=str(theField),fixuvw=False, direction=skyPosAtField[theField]) # Change flags in FIELD table to signalize that directions has been set # - the values in the FIELD table are not relative offsets anymore # but the absolute RA/DEC coordinates in ICRS (ref. frame id=21 <=> ICRS) # (NB: right after fixplanets() call it was id=0 <=> J2000). tb.open(msc+'/FIELD', nomodify=False) for sf in SunFields: tb.putcell("RefDir_Ref", sf, 21) tb.putcell("DelayDir_Ref", sf, 21) tb.putcell("PhaseDir_Ref", sf, 21) tb.close() #========================================================================= print("\n\n******* Calibration finished.*******\n\n") #===========================================================================