#! /usr/bin/env python # -*- coding: iso-8859-1 -*-# # # ************** Solar ALMA Data Reduction Script Generator ************** # # Single-purpose code for generation of solar-specific calibration # and imaging scripts. To be used by QA2 analysts at ARCs. # The idea behind is to have a similar generator as the one by Eric Villard for # non-solar data. Current implementation is, however, much less robust (right # ASDM data are expected to be present), and, as written, less general (just # solar science and CSV/EoC ALMA test INT data from 7m+12m array). # # begin: 2019-03-14 # last update: 2025-12-12 # Miro Barta , EU ARC - Czech node # # How to use: # # 1. import SolSG as ssg # 2. sso = ssg.stuffForSolarDataReduction() # 3a. sso.generateSolCalScript(asdmName,doQA2=True) # 3b. sso.generateSolImgScript(asdmName,doQA2=True) #----------------------------------------------------------------------------- #--- Import general libs import os import numpy as np #--- Import selected CASA tools and tasks. from casatools import msmetadata as msmdtool from casatools import quanta as qatool from casatools import table as tbtool from casatools import measures as metool from casatasks import importasdm #--- Import analysis utils import analysisUtils as aU es=aU.stuffForScienceDataReduction() #============================================================================= # class scienceTarget #============================================================================= class scienceTarget: """ Attributes & methods for the science target (the Sun). Attributes: allFields (list of ints): List of all target fields. May contain purelly OFF_SOURCE field 0 that is attached to ephemeris and serves as a reference to which the science fields in a mosaic have offsets in the FIELD table. sciFields (list of ints): Just the science fields (the intent OBSERVE_TRAGET#ON_SOURCE) refField (int): The reference field attached to ephemeris. refTime (CASA qanta): The time when the first science field started to be scanned. For exactly this time the ephemeris is interpolated and RA/DEC found by this interpolation are set to the reference field. It is also the same date/time that is later written as DATE_OBS keyword in the header of final product FITS image - can be used for transformation to helioprojective X,Y. fieldsDirections (dictionary of measures): Ra/Dec attributed to each field at the reference time. Methods: bool isMosaic(): Returns T if the target is a mosaic, F if single-pointing. string targetFields(): Returns the all target fields as a CASA-string. void printInfo(): Prints information about science target to terminal. """ def __init__(self,obsMsmd): """ Constructor. @param obsMsmd: MS meta-data. """ #--- Explore the MS metadata, extract the science-target info. allFieldsArray=obsMsmd.fieldsforintent("*OBSERVE_TARGET*") sciFieldsArray=obsMsmd.fieldsforintent("*OBSERVE_TARGET#ON_SOURCE*") sciScan=obsMsmd.scansforintent('*OBSERVE_TARGET#ON_SOURCE*') sciScanTimes=obsMsmd.timesforscan(sciScan[0]) #--- Make extracted metadata to be the science-target attributes. self.allFields=allFieldsArray.tolist() # Python list self.sciFields=sciFieldsArray.tolist() # Python List self.refTime=qatool().quantity(sciScanTimes[0],'s') # qa (value & unit) #--- Find the phase-center (reference) field using the FIELD table # Open table with phase-center offsets. fTable=tbtool() fTable.open(obsMsmd.name()+'/FIELD', nomodify=True) phsCenOff = fTable.getcol("PHASE_DIR") fTable.close() # RA/DEC offsets for all the solar fields SunRA_offs=phsCenOff[0][0][self.allFields] SunDEC_offs=phsCenOff[1][0][self.allFields] # Find the field with minimum offset as a candidate dist2CenSqr=SunRA_offs**2+SunDEC_offs**2 centerIndex=np.argmin(dist2CenSqr) fieldRef=self.allFields[centerIndex] # The reference field should have zero PHASE_DIR offset in the # FIELD table before the fixplanets() is being called. # Check, whether it is true. Should this test fail, # the analyst must set the 'fieldRef' manualy # This is not 100% 'standard': Warn and let the analyst to decide. if(phsCenOff[0][0][fieldRef]!=0.0 or \ phsCenOff[1][0][fieldRef]!=0.0): print("Reference field candidate: "+str(fieldRef)+".\n") print("Candidate for the reference field has non-zero offset:\n") print("RA offset = "+str(phsCenOff[0][0][fieldRef])) print("DEC offset = "+str(phsCenOff[1][0][fieldRef])) print("Possible issue with setting absolute Ra/Dec coordinates.") print("Check FIELD table and try to find manual workaround.\n") # Set the found reference-field candidate as the sci-target attribute. self.refField=fieldRef #--- Find the Ra/Dec coordinates of all solar field at the # reference time via msmd.phasecenter() # msmd.phasecenter() needs the reference time expressed as a MJD and # specified as a measure::epoch object, including time-frame (UTC). refTimeMJD=metool().epoch('utc',self.refTime) # Initialise an empty map (dictionary) self.fieldsDirections=dict() # Fill it field-by-field. for field in self.allFields: # Find the sky position (RA/Dec) of a given field at ref. time... fDir=obsMsmd.phasecenter(field,refTimeMJD) # ... and append it to the attribute holding a list of all positions self.fieldsDirections[field]=fDir #------------------------------------------------------------------------- def isMosaic(self): """ Returns True if the target is a mosaic, False if single-pointing. """ nSciFields=len(self.sciFields) return (nSciFields >1) #------------------------------------------------------------------------- def targetFields(self): """ Returns a string to be passed as a 'field' argument in applycal() and tclean(). Note: We have to include the field 0 (in case of mosaic) whether it is associated with a scicence scan or not. The 'phantom' field=0, albeit may not even be present in the mosaic as a science target field (MS.split.cal contains no data associated with this field), it is still in the FIELD table and THIS is the field that is clipped to ephemeris (even if it is purelly OBSERVE_TARGET#OFF_SOURCE intent!) It is usually located in the mosaic center, but need not be. """ if (self.isMosaic()): nSunFields=len(self.allFields) tgtFields="0, "+str(self.allFields[1])+"~" \ +str(self.allFields[nSunFields-1]) else: tgtFields='0' return tgtFields #------------------------------------------------------------------------- def printInfo(self): """ Print info/summary about the target for the user. """ if self.isMosaic(): print("Your observation is a mosaic.") # How many science fields do we have in the mosaic? nSciFields=len(self.sciFields) print("Number of mosaic pointings is "+str(nSciFields)) # More details. print("First mosaic field is " + str(self.sciFields[0])) print("Last mosaic field is " + str(self.sciFields[nSciFields-1])) print('All target fields = '+'\"'+self.targetFields()+'\"') else: print("Your observation is a single-pointing.") print("The science target field id is 0.\n") # Create a quanta tool for conversions. convertor=qatool() # Convert/export the refTime quantity to a string and print it. refTimeStr=convertor.time(self.refTime,form="ymd",prec=8)[0] print('Reference time: '+refTimeStr) # Print the reference field. print('Reference field: '+str(self.refField)) # Print the reference (phase center) Ra/Dec exported as a string. RA_Cen=(self.fieldsDirections[self.refField])['m0'] # QA object DEC_Cen=(self.fieldsDirections[self.refField])['m1'] # QA object # RA/DEC as a formated string in HMS/DMS format RA_Cen_hms=convertor.angle(RA_Cen,prec=10,form='time')[0] DEC_Cen_dms=convertor.angle(DEC_Cen,prec=10)[0] print('Reference field (phase center) RA/Dec = '+RA_Cen_hms \ +' / '+DEC_Cen_dms+"\n") # Close the convertor. convertor.done() #============================================================================= # class calibrators #============================================================================= class calibrators: """ Encapsulation of calibrators. Attributes: fluxCalField (string): Field for intent CALIBRATE_FLUX. phaseCalField (string): Field for intent CALIBRATE_PHASE. bpassCalField (string): Field for intent CALIBRATE_BANDPASS. ampCalField (string): Field to be actually used in setjy() - may be bpassCalField or fluxCalField. bpassCalScan (string): Which scan has been used for B-pass calibration. ampCalSrcName (string): Source-name for actual flux calibrator. """ def __init__(self,obsMsmd): """ Initialize the calibrator attributes """ #--- Find the fields of specific calibrators (by intent). # Flux calibrator: This one need not to be present. # And even if the object with the intent FLUX exists, it may be not # the right calibrator - see self.fluxCalQualityMetrics(). fcFields=obsMsmd.fieldsforintent("*FLUX*") if(len(fcFields)>0): self.fluxCalField=str(fcFields[0]) else: self.fluxCalField='' # Phase calibrator: always present self.phaseCalField=str((obsMsmd.fieldsforintent("*PHASE*"))[0]) # Band-pass calibrator: always present self.bpassCalField=str((obsMsmd.fieldsforintent("*BANDPASS*"))[0]) self.bpassCalScan=str((obsMsmd.scansforintent("*BANDPASS*"))[0]) #--- Which calibrator we shall use for flux-density in setjy()? # NB: So far we select the Band-Pass calibrator but in future it # should be checked, whether the designated flux calibrator - if # present - is better than the band-pass one via (insofar not yet # implemented) self.fluxCalQualityMetrics(). # if(self.fluxCalField==''): self.ampCalField=self.bpassCalField else: """ This comes into play after self.fluxCalQualityMetrics() gets implemented! --- qFluxCal=self.fluxCalQualityMetrics(self.fluxCalField) qBpCal=self.fluxCalQualityMetrics(self.bpassCalField) if(qFluxCal > qBpCal): self.ampCalField=self.fluxCalField else: self.ampCalField=self.bpassCalField --- For now use the B-pass calibrator for setjy() as default and let the user decide to set manually the (designated) F-calibrator. """ self.ampCalField=self.bpassCalField # This is not 100% 'standard': Warn and let the analyst to decide. print('WARNING: Designated flux calibrator found at field '\ +self.fluxCalField+".\n") print('Check its flux and last measurement using aU.getALMAFLux()') print('Use it for setjy() and subsequent fluxscsale() only if it ') print("is better than the (default) band-pass calibrator!\n\n") #--- We need to find selected amp-cal source name for aU.getALMAFlux() # There is no direct msmd method for that - we have to work around # and search in the SOURCE table. # Get all (non-unique) source IDs (a long numpy.ndarray) allSrcIds=obsMsmd.sourceidsfromsourcetable() # Find the source ID for the band-pass calibrator. ampCalSrcID=obsMsmd.sourceidforfield(int(self.ampCalField)) # Find the 1st index of the amp-cal in the array of all src IDs. # This means the 1st row in SOURCES table where the amp-cal occurs. ampCalRow=list(allSrcIds).index(ampCalSrcID) # Now read the source name from the source table. self.ampCalSrcName=(obsMsmd.sourcenames())[ampCalRow] #------------------------------------------------------------------------- def allCalsStr(self): """ All calibrators (a unique set) as a single comma-separated string. To be used in gaincal(field=allCals) """ allCalsString=self.bpassCalField if(self.fluxCalField!='' and self.fluxCalField!=self.bpassCalField): allCalsString+=','+self.fluxCalField if(self.phaseCalField!=self.bpassCalField and self.phaseCalField!=self.fluxCalField): allCalsString+=','+self.phaseCalField return allCalsString #------------------------------------------------------------------------- def allCalsLst(self): """ All calibrators (a unique set) as a Python list (of strings). To be used in a loop over applycal(field=allCals[i]). """ allCalsList=[self.bpassCalField] if(self.fluxCalField!='' and self.fluxCalField!=self.bpassCalField): allCalsList+=[self.fluxCalField] if(not(self.phaseCalField in allCalsList)): allCalsList+=[self.phaseCalField] return allCalsList #------------------------------------------------------------------------- def printInfo(self): """ Print info on calibrators. """ print('Bandpass & flux calibrator found at field '+self.bpassCalField) print("Flux density calibrator source name: "+self.ampCalSrcName) print('Phase calibrator found at field: '+self.phaseCalField+"\n") #------------------------------------------------------------------------- def fluxCalQualityMetrics(fCalField): """ Returns a single float number as a measure of quality for the flux calibrator. The quality is based on (1) flux density and, (2) time separation between the last measurements in the ALMA calibrator database and the date of observations. Rationale: Quite often we have declared (primary) flux calibrator (with intent CALIBRATE_FLUX) in solar science data which, however, is either weak or its measurement in database is old. Hence, the band-pass calibrator is frequently use as a primary calibrator: The function sol_calflux2() even assumes this as a standard. This method should calculate a score (the higher, the better) to help to decide whcih calibrator should be used for setjy() and subsequent fluxscale() """ # TO DO return 0.0 #============================================================================= # class observation #============================================================================= class observation: """ Basic minimum data about the observation that are necessary for data reduction & imaging. Contains the children objects scienceTarget and calibrators plus a few further service data. Attributes: theSun (class scienceTarget): Encapsulated properties of the target theCalibrators (class calibrators): Encapsulate properties of calibrators. band (int): Observing band. refFreq (CASA quanta): Representing (reference) observing frequency. obsDate (string): Date of obsrevation in YYYY-MM-DD format. """ def __init__(self,asdmName): #--- Import the ASDM msName=asdmName+'.ms' if not os.path.exists(msName): # MS has not been imported yet, do it now. print("Importing ASDM to MS...") importasdm(asdm = asdmName, vis = msName, lazy = True, asis='Antenna Station Receiver Source CalAtmosphere ' +'CalWVR CorrelatorMode SBSummary CalDevice') # Check whether the MS is here now (it should be). if not os.path.exists(msName): # ASDM import failed for some reason - exit now. print("ASDM import failed, exiting!\n\n") return #--- The MS exists, let us explore it via its meta-data tool. print("Exploring basic MS parameters...") obsMsmd=msmdtool() obsMsmd.open(msName) #--- Specific attributes: Science target and calibrators self.theSun=scienceTarget(obsMsmd) self.theCalibrators=calibrators(obsMsmd) #--- Further general attributes at which the data reduction depends: # In which Band and at which reference frequency are we observing? # What is the observation date? # The science SPWs for solar data are 5,7,9, and 11. # Should there be a change, one can use # scispws=msmd.almaspws(tdm=True,fdm=True) in conjunction with # msmd.spwsforintent("OBSERVE_TARGET#ON_SOURCE") - the intersection # of the two lists are the science SPWs. # Let us be simple for now and use SPW 7 for ALMA Band identification. refSPW=7 # The aU.bandforspw(spw,mymsmd) uses msmd.namesforspws() internally. self.band=aU.bandforspw(refSPW,mymsmd=obsMsmd) # Find the reference observing frequency (in the reference SPW). self.refFreq=(obsMsmd.reffreq(refSPW))["m0"] # As quantity # Observing date to be extracted form the reference time for the Sun. self.obsDate=qatool().time(self.theSun.refTime, form=["fits","no_time"])[0] # MS exploration finished. obsMsmd.done() #------------------------------------------------------------------------- def rfStr(self): """ Returns the reference frequency as a string (value + unit). """ return str(self.refFreq["value"])+str(self.refFreq["unit"]) #------------------------------------------------------------------------- def printInfo(self): """ Prints basic info about the observation. """ print("\nYour observation is in Band "+str(self.band)+".") print("The reference frequency is "+self.rfStr()+".") self.theSun.printInfo() self.theCalibrators.printInfo() #============================================================================= # class stuffForSolarDataReduction #============================================================================= class stuffForSolarDataReduction: """ Encapsulation for Solar ALMA Data Reduction Utils""" # Note: As a part of good object-programming practice it would be nice to # encapsulate the routines in sun_redcution_util.py as static methods # of this class, too. So far they are imported from an external module #------------------------------------------------------------------------- def generateSolCalScript(self,asdmName,doQA2=True): """ Generates the script for solar INT data calibration (single EB) Implementation is very simple: It reads the static template (the 'sceleton', which remains constant for all solar INT EBs) and adds/replaces the lines with the EB-specific code. Result is written into the 'standard' solar calibration script. Inputs: asdmName (string): Name of the ASDM to reduce doQA2 (bool): Generate the script for the QA2 analyst (True) or for the archive (False). """ #--- Import MOUS, extract observation info required for data reduction. print("\nImporting and exploring the MOUS.") myObs=observation(asdmName) # Print a brief observation summary for the analyst myObs.printInfo() #********************************************************************** #--- Generate the code lines (snippets) that shall replace the KEYs # strings in the code template. Mapping between the replacable keys # in the script template and real lines of code in the final # generated script is provided by a dictionary val=codeSnipets[key] # Generated text is a Python code, which requires proper indentation. indL1=' ' # Indent level 1 indL2=' ' # Indent level 2 indL3=' ' # Indent level 3 # Create an empty dictionary ... codeSnipets=dict() # ... and append the key/value (code replacements) pairs to it: #--------------------------------------------------------------------- # A piece of code to replace the __KEY_ASDM_NAME__ line in template ... codeSnipets['ASDM_NAME']="asdm = \'"+asdmName+'\'' #--------------------------------------------------------------------- # ... and so on with the other lines to be replaced. # Define reference antenna. # Autodetection via aU.getRefAntenna() may fail - protect it by # a try-catch block. try: refAnt='\''+es.getRefAntenna(asdmName+'.ms')+'\'' except IndexError: print('!!! Auto-detect reference antenna failed !!!') print('The string \'refAnt\' in the calibration script set empty.') print('Has to be reset by the analyst latest at Step 13 !!!\n') refAnt='\' \' # MUST be set manually!' codeSnipets['REF_ANT']="refAnt="+refAnt #--------------------------------------------------------------------- # Assign fields/scans for calibrators & science target codeSnipets['OBS_FIELDS']="fluxCalField = \'" \ +myObs.theCalibrators.fluxCalField+'\'\n' \ +"phaseCalField = \'" \ +myObs.theCalibrators.phaseCalField+'\'\n' \ +"ampCalField = \'" \ +myObs.theCalibrators.ampCalField+'\'\n' \ +"bpassCalField = \'" \ +myObs.theCalibrators.bpassCalField+'\'\n' \ +"allCals = \'" \ +myObs.theCalibrators.allCalsStr() \ +"\' # All calibrators (fields) as a comma-separated string.\n" \ +"allCalsList = " \ +str(myObs.theCalibrators.allCalsLst()) \ +' # All calibrators (fields) as a Python list. \n' \ +"bpassCalScan = \'" \ +myObs.theCalibrators.bpassCalScan \ +"\' # Id of the scan for the bandpass calibration.\n" \ +"SunFields = " \ +str(myObs.theSun.allFields) \ +' # All Sun fields as a list (of strings).\n' \ +"targetFields = \'" \ +myObs.theSun.targetFields() \ +"\' # All Sun fields as a string.\n" #--------------------------------------------------------------------- # Extra diagnostic plot of listobs in the QA2 mode. codeSnipets['LISTOBS_PLOT']=indL1+"# Observing summary as a plot.\n" \ +indL1+"es.listobs3(mso, figfile=asdm+'_listobs3.png')" if(doQA2) \ else "" #--------------------------------------------------------------------- # Extra diagnostic plots for Tsys in the QA2 mode. # !!!NB: Replace plotcal() by CASA6-compatible command tSysDiagQA2=indL1+"# Make the plots for Tsys checks in the QA2 mode.\n" tSysDiagQA2+=indL1+"es.checkCalTable(mso+'.tsys', msName=mso, " \ +"interactive=False)\n" \ +indL1+"ant_names = getAntennaNames(mso + '.tsys')\n" \ +indL1+"os.system('rm -rf ' + mso + '.tsys' + '.plots_ext')\n" \ +indL1+"os.system('mkdir ' + mso + '.tsys' + '.plots_ext')\n\n" \ +indL1+"for ant_nm in ant_names:\n" \ +indL2+"plotms(mso+'.tsys',xaxis='time',yaxis='tsys',\n" \ +indL3+"antenna=ant_nm,gridrows=2,gridcols=2,iteraxis='spw',\n" \ +indL3+"plotfile=mso+'.tsys'+'.plots_ext/timeprof_'+ant_nm+'.png'" \ +",\n"+indL3+"coloraxis='corr',overwrite=True,showgui=False)\n\n" \ +indL2+"plotms(mso+'.tsys',xaxis='chan',yaxis='tsys',\n" \ +indL3+"antenna=ant_nm,gridrows=2,gridcols=2,iteraxis='spw', \n" \ +indL3+"plotfile=mso+'.tsys'+'.plots_ext/specprof_'+ant_nm+'.png'" \ +",\n"+indL3+"coloraxis='corr',overwrite=True,showgui=False)\n\n" codeSnipets['TSYS_PLOT']=tSysDiagQA2 if(doQA2) else "" #--------------------------------------------------------------------- # Generate the antenna-position correction table # This will generate the string with gencal(caltype = 'antpos',...) # preceded by removal of the old caltable. apcCommand=es.correctMyAntennaPositions(asdmName+'.ms',verbose=False, checkForAntennaMoves=False) # Insert the gencal() command generated above only if it is non-empty. if(apcCommand): # Get rid of the OS command removing the old antPos caltable - we # have it already in the cal-script template, start with gencal(). genCalBeg=apcCommand.index('gencal(') genCalCommand=apcCommand[genCalBeg:] # Replace actual MS name with its shortcut defined in the script. # Have to do it in two steps: genCalCommand=genCalCommand.replace('\''+asdmName+'.ms.antpos\'', 'mso+\'.antpos\'') genCalCommand=genCalCommand.replace('\''+asdmName+'.ms\'','mso') else: # No gencal() command has been generated. This may be due to old # AIV/science/PadData database but it may happen even with the # most-recent one. We just put a comment into the script. genCalCommand='# WARNING: No gencal() command generated by ' \ +'es.correctMyAntennaPositions().\n' \ +'# Check the SSG on-screen print-outs or the CASA log.\n' # Insert the command into the dictionary with a proper indentation. codeSnipets['ANT_POS_TABLE']=indL1 \ +genCalCommand.replace('\n','\n'+indL1) #--------------------------------------------------------------------- # Applycal for Tsys and TantPos to all calibrators # Note: We are using gainfield='nearest' in order to handle all # cal-tors homogeneously, either they have their own Tsys scans # or not. applyPreCalibs='' for caltor in myObs.theCalibrators.allCalsLst(): applyPreCalibs+=indL1+"applycal(vis = mso,\n" \ +indL3+"field = \'"+str(caltor)+"\',\n" \ +indL3+"spw = \'5,7,9,11\',\n" \ +indL3+"gaintable = [mso + \'.tsys\', mso+\'.antpos\'],\n" \ +indL3+"gainfield = [\'nearest\',''],\n" \ +indL3+"interp = [\'linear,linear\', \'nearest\'],\n" \ +indL3+"spwmap = [tsysmap,[0]],\n" \ +indL3+"calwt = [True, False],\n" \ +indL3+"flagbackup = False)\n\n" codeSnipets['APPLY_PRECALIBS']=applyPreCalibs #--------------------------------------------------------------------- # Plot the calibration weights statistics in the QA2 mode. codeSnipets['CAL_TANT_PLOT']=indL1 \ +"# Plot the solar-scan antenna temperatures Tant.\n" \ +indL1+"sol_tsunplot(asdm)\n\n" \ +indL1+"# Plot the cal. weights statistics.\n" \ +indL1+"es.getCalWeightStats(mso)" if(doQA2) else "" #--------------------------------------------------------------------- # Scans with which intent shall continue for further processing after # splitting-out just the science data? sciIntents='*BANDPASS*,*PHASE*,OBSERVE_TARGET#ON_SOURCE' \ if(myObs.theCalibrators.fluxCalField=='') else \ '*FLUX*,*BANDPASS*,*PHASE*,OBSERVE_TARGET#ON_SOURCE' codeSnipets['SPLIT_INTENT']=indL1+'sciScansIntent = \''+sciIntents+'\'' #--------------------------------------------------------------------- # Flagging of channels contaminated by the Earth atmospheric lines flagAtmLines=indL1 if(myObs.band == 6): flagAtmLines+="# BB2: 231.3GHz\n" \ +indL1+"# BB4: 247.7GHz (weak), 248.2 (strong)\n\n" flagAtmLines+=indL1+"# BB2\n" \ +indL1+"flagdata(vis = mss,\n" \ +indL3+"spw=\'1:105~114\',\n" \ +indL3+"mode = \'manual\',\n" \ +indL3+"flagbackup = False)\n\n" flagAtmLines+=indL1+"# BB4 - weak line \n" \ +indL1+"flagdata(vis = mss,\n" \ +indL3+"spw=\'3:40~55\',\n" \ +indL3+"mode = \'manual\',\n" \ +indL3+"flagbackup = False)\n\n" flagAtmLines+=indL1+"# BB4 - strong line \n" \ +indL1+"flagdata(vis = mss,\n" \ +indL3+"spw=\'3:68~84\',\n" \ +indL3+"mode = \'manual\',\n" \ +indL3+"flagbackup = False)\n\n"\ elif(myObs.band == 7): flagAtmLines+=indL1+"# BB3\n" \ +indL1+"flagdata(vis = mss,\n" \ +indL3+"spw='2:90~126',\n" \ +indL3+"mode = 'manual',\n" \ +indL3+"flagbackup = False)\n\n" flagAtmLines+=indL1+"# BB4\n" \ +indL1+"flagdata(vis = mss,\n" \ +indL3+"spw='3:0~30',\n" \ +indL3+"mode = 'manual',\n" \ +indL3+"flagbackup = False)\n\n" else: # !! Check that we really do not have ATM lines in Band 5 !! flagAtmLines+="# No contaminating ATM lines in Band " \ +str(myObs.band)+" to flag." codeSnipets['FLAG_ATM_LINES']=flagAtmLines #--------------------------------------------------------------------- # Flux calibration: Set Jy to the flux calibrator # Query the catalog for the calibrator properties. fcVals=aU.getALMAFlux(sourcename=myObs.theCalibrators.ampCalSrcName, date=myObs.obsDate, frequency=myObs.rfStr(),silent=True) # Extract just the flux-density and spectral slope from all the values. fluxDensity=str(fcVals["fluxDensity"]) spSlope=str(fcVals["spectralIndex"]) # Compose the setjy() command. fluxCalCMD=indL1+"setjy(vis = mss, field = ampCalField, \n" \ +indL3+"spw = \'0,1,2,3\',\n" \ +indL3+"standard = \'manual\',\n" \ +indL3+"fluxdensity = ["+fluxDensity+" ,0, 0, 0],\n" \ +indL3+"spix = "+spSlope+",\n" \ +indL3+"reffreq = \'"+myObs.rfStr()+"\')" # Define the flux-density calibration command. codeSnipets['SET_JY']=fluxCalCMD #--------------------------------------------------------------------- # Diagnostic plots for band-pass calibration (Step 13) in the QA2 mode. bpCalDiag=indL1 \ +"# Plot the pre-bandpass and bandpass calibration tables.\n" bpCalDiag+=indL1+"es.checkCalTable(mss+'.ap_pre_bandpass', " \ +"msName=mss, interactive=False)\n" bpCalDiag+=indL1+"es.checkCalTable(mss+'.bandpass', " \ +"msName=mss, interactive=False)\n\n" codeSnipets['BP_CAL_PLOT']=bpCalDiag if (doQA2) else "" #--------------------------------------------------------------------- # Diagnostic plots for gain calibration (Step 15) in the QA2 mode. gainCalDiag=indL1+"# Plot the calibration tables of:\n" \ +indL1+"# - Itegration-based phases, \n" \ +indL1+"# - scan-based phases, \n" \ +indL1+"# - scan-based amplitudes (rel.), and\n" \ +indL1+"# - scan-based fluxes (flux-scaled).\n\n" gainCalDiag+=indL1+"es.checkCalTable(mss + '.phase_int', " \ +"msName=mss, interactive=False)\n" gainCalDiag+=indL1+"es.checkCalTable(mss+'.phase_inf', "\ +"msName=mss, interactive=False)\n" gainCalDiag+=indL1+"es.checkCalTable(mss + '.ampli_inf', "\ +"msName=mss, interactive=False)\n\n" gainCalDiag+=indL1+"es.fluxscale2(caltable = mss+'.ampli_inf',\n" \ +indL3+"removeOutliers=True,\n" \ +indL3+"msName=mss, writeToFile=True,\n" \ +indL3+"preavg=10000)\n\n" # Insert the gain-calibration diagnostic into the code in the QA2 mode. codeSnipets['G_CAL_PLOT']=gainCalDiag if (doQA2) else "" #--------------------------------------------------------------------- # Sky positions (RA/Dec) for the target fields (Step 20) # We need to convert them from a list of measures to the list of strings # Initialise an ampty list (of strings). sunFieldDirs=[] for theField in myObs.theSun.allFields: fieldPosition=myObs.theSun.fieldsDirections[theField] # RA/Dec in radians (quanta: Value + unit) RA=fieldPosition['m0'] DEC=fieldPosition['m1'] # Convert the values in radians to usual HMS/DMS scale. Prepend the # epoch/reference frame J2000 to the string. fPosStr='J2000 '+aU.rad2radec(RA['value'],DEC['value'], verbose=False, prec=1, hmsdms=True, delimiter=' ') sunFieldDirs.append(fPosStr) # Make the list to be a string and improve its readability. sfdStr=str(sunFieldDirs).replace(', ',',\n ') # Insert a paparmeter 'skyPosAtField' with a value sfdStr to the code. codeSnipets['SUN_FIELD_DIRS']=indL1+"skyPosAtField="+sfdStr+"\n" # Dictionary completed #********************************************************************** # Compose & write the calibration script. print('Writing the calibration script...') #--- Open a file with the static script template solCalPath="/home/barta/ALMA/QA2/science/qa2/solar/" scrTemplate=open(solCalPath+"SSG/solCalScriptTemplate_INT.py","r") #--- Parse the script template/skeleton, replace the __KEY_xxx__ # with a MS-taylored lines of code # Calibration script code to be stored in this string line-by-line. # Initialise it as an empty string for now. calScriptCode="" # Read the template line-by-line, search for the lines with __KEY_xxx__ for line in scrTemplate: if "__KEY_" in line: # Split the line to continguous stripped parts separated by # white-spaces. The full-key (including the 'KEY' identifier # and all uderscores) is at the 1st position. fullKeyId=line.split()[0] # Strip the prefix '__KEY_' and suffix '__'. theKey=fullKeyId[6:-2] # Replace the key with the specific snippet of the code. calScriptCode+=codeSnipets[theKey]+'\n' else: # This is a "static" part of the code, just copy the template. calScriptCode+=line # Template parsing finished. scrTemplate.close() #--- Open the file with a script, write the contents & close scriptFileName=asdmName+'.ms.scriptForCalibration.py' scriptFile=open(scriptFileName,'w') bytesWritten=scriptFile.write(calScriptCode) scriptFile.close() print("\nScript for calibration written to "+scriptFileName \ +" ("+str(bytesWritten)+" bytes).\n") #-------------------------------------------------------------------------