Useful information for navigating all these scripts:

• fluxes.py and run_gainCheck.csh are in /array/obs/bin/
• quality is in /opt/rt/scripts/

Updates:

• Amber dropped in 11/22 version of fluxes_temp.py on Friday, 11/25 coinciding with her fluxcal run - DONE on 11/25 at 12:15pm - catalog files in fluxes_temp.py incremented to v2 and old script catalog log file updated to oldstyle3

Logging new vs. old script outputs:

• new and old scripts write out to log files in ~/chat/fluxcal/ on high site
• Updates:
  • Oct 21: Amber changed new script to write to sci1 and sci2 logfiles separately as we were getting incorrect logging due to sci1/sci2 conflicts in writing to the log file; incremented log files
  • Nov 25: Amber made major edits to the new script; incremented log files
  • Nov 29: Amber fixed a bug with config_wb in processFile(); updated on high site

• Make sure running well at Illinois
• Better planet models?

1. Extract spectral channels only (from original vis file) → speconly.mir
2. Remove bad windows for sci2/CARMA23 obs (from speconly.mir) → vis_sci2.mir
3. Just wide-band windows, pol(LL) (from speconly.mir or vis_sci2.mir) → wb.mir
4. Amplitude flagging creates intermediate files described below, but flags wb.mir. (These steps made necessary because smauvplt cannot average if there is a planet in the observation…reason unknown)
   1. planet-only uvmodel options=divide → wb_planet.mir
   2. uvaver wb_planet.mir → wb_planet_avg.mir
   3. calibrators w/out planet → wb_noplanet.mir
5. Passband calibration → pb.mir
6. Phase selfcal planet only → primary.mir
7. Phase selfcal each calibrator → source.mir

Currently, two log files are maintained - fluxcal.log and fluxcal.cat. Fluxcal.cat is the one that is to be put into FluxSource_newadd.cat (currently though, the file names are written out at the start of each line, for testing purposes). Fluxcal.log records what happened when a track does or does not successfully complete. Failure modes include:

• notScience: track does not start with c0, c1 or c2
• Len/Grade: track is not A or B, length > 1
• NoPrimCal: track does not contain a primary calibrator (Mars, Uranus or Neptune)
• noConfig: track contains 0 or > 1 correlator configurations
• noWBwin: track does not contain any wide band windows
• allAmpFlag: too much bad data - all of it was flagged in the amplitude flagging section
• configProb0/1/2: number of corr configs != 1 in original file / after wb win extraction / after calibrator extraction (first two in flux, last in processFile)
• passband: unable to get successful passband calibration for any refant (success here = refant is good and at least 3 good antennas)
• noGoodFluxes: after loop through cals running bootflux, unable to extract good fluxes for some reason
• planetRes: bootflux crashed with a fatal error, likely because the planet is resolved
• bootfluxVar: the bootflux values varied by more than 25% for that calibrator (in this case, each offending calibrator will be output with this error to fluxcal.log with the number of time point indicated in parenthesis next to the calibrator name)

def flux(vis=None, db=False, plotResult=False, passtrack=False, ampFlag=2)

This is the main driving program: this is what you call to calculate fluxes for a track, ie

import fluxes_temp 
fluxes_temp.flux('blah.mir',$subarray)

Currently, this is called from run_gainCheck.csh in /array/obs/bin/ at the high site.

• 11/29 AB added second call to run 1cm PACS dataset if the observation is PACS since PACS observations are run in 'sciall' (similar to sci1), so that without this special call AB added, the 1cm PACS datasets will not be run through fluxes_temp since there are no sci2 observations.

Steps of analysis:

1. Setup:
   • check vis file exists
   • check that dataset is a science track - must start with cx, c0, c1, c2
   • run listobs (to logfile listobs.log) and calculate grade, track length, antennasOnline and check for primary calibrators
2. check that track has primary cal - only move on if contains MARS, NEPTUNE or URANUS (also, checkForPrim checks that there is a maximum of 1 planet in the observation)
3. check grade/tracklength conditions: move on if grade=A or B and tracklen > 1
4. figure out obsmode from antennasOnline list: CARMA15, CARMA23 or sci2
5. due to issue with wideband channels being messed up in some sci2 tracks, first extract only the spectral channels with uvcat options=nowide, writing speconly.mir.
6. get correlator configuration of original vis file with getCorrelatorConfigurations
   CRITERIA: there must be at 1 wide-band window - if 0, failure mode noConfig; if > 1, failure mode configProb0
7. remove bad windows from 3.5m data (sci2 or CARMA23 tracks): if 1cm (sci2 only), remove windows 17-32, if 3mm, remove windows 1-16 for sci2, windows 1-8 for CARMA23. Write vis_sci2.mir. Get updated correlator configuration.
8. get nchan and min/max freq of correlator config, make list of wideband windows
9. create file of just wideband data (from calibrators.mir) selecting windows and min,max freq and only LL polarization. Writes wb.mir
10. read new correlator configurations with getCorrelatorConfigurations
11. write visfile name to temporary recording file (chat/fluxcal/fluxsource_temp_sci1/2.cat)
12. loop through bands (1cm, 3mm, 1mm):
    • check if any of the configurations lie within the band
    • if so, process the data if the number of configurations is 1:
      processFile(calibrators.mir,configurations[nconfig],sci2,obsblock=obsblock, db=db)
    • if no configurations in this band, output 'No data in …'
    • if multiple configurations, failure mode configProb1: output 'Not processing vis= for the b band: multiple correlator configurations'

def calibratePassband(mirfile, calibrators, refant, antennasOnline)

1. standard mfcal (with all cals) with INTERVAL_PB (currently set to 0.5) and gpedit options=phase sets amp gains to 1:

       mfcal vis=mirfile interval=INTERVAL_PB refant=refant
       gpedit vis=mirfile options=phase

2. calculate signal to noise on each calibrator, finding max snr
3. passband calibrate with max snr calibrator, same mfcal line as before but with select= and no gpedit line; checks logfile for FatalError with containsFatalError function
4. print median and rms gains from gplist command
   write gain median and rms files to dict's called gains and rms for those antennas which are in antennasOnline array
   if gain == 0 (or conversion of gain value to float failed), add antenna to badAntennas array
5. determine good ranges for antenna gain and rms values, appending those w/in good range to goodAntennas list (others append to badAntennas)
   goodrange for gains = median(median gains)*(1 +- 0.3); for rms = max of 2*median(rms vals) or 0.1
6. success is True if refant is in goodAntennas list AND there are at least 3 goodAntennas. If successful, writes out pb.mir file (removing badAntennas):

       rm -rf pb.mir
       gpedit vis=mirfile options=phase
       uvcat vis=mirfile out=pb.mir options=nocal select=-badAntennas

7. returns success,mirout = True,pb.mir or False,''

def flagMiriadFile(mirfile, correlator, allcals, ampFlag=2)
mirfile = name of miriad file
correlator = a dictionary of dictionaries describing correlator windows (ie. 'nchan')
allcals = list of calibrators to flag
ampFlag = 0,1 or 2 (default = 2) for no amplitude flagging, just relative amplitude flagging (based on median amp) or both relative and absolute amp flagging

Applies flagging to mirfile for the following things:

1. shadowing: csflag vis=
2. high elevation - uvflag select=el(85,90)
3. edge channels - first determine min/max number of channels in correlator windows specified in correlator variable
   • if all windows are same nchan (and !=15), just do uvflag vis= flagval=f edge=EDGECHAN_FLAG (3)
   • if not, flags each window individually
   • comment here says that when uvflag bug is fixed, this can be simplified - ask Peter again if this has been fixed!
4. spuriously large amplitudes:
   • if there is a planet, uvcat to take out planet → wb_noplanet.mir
   • uvindex on wb.mir → wb_uvindex.log
   • go through each source, keeping track if a window is EVER marked bad:
     • planet: ampFlagPlanet in each window
     • non-planet: calcAmpFlagging in each window, comparing avg window fluxes at the end - any windows outside +/- 40% from median are marked bad
   • flag bad windows (bad in any source → flag for entire obs)

def processFile(mir_wideband, config_wb, obsmode, obsblock=None, db=False, plotResult=False, ampFlag=2)

1. get source list and calibrators list (including planets)
2. flagMiriadFile(wb.mir or wb_sci2.mir, config_wb, calibrators list)
3. extract only calibrator data → calibrators.mir
4. get updated correlator configuration (config_wb) and calculate mean freq of windows (freq)
5. get potential passband calibrator list (selecSources planets=False)
6. check for primary calibrators (PRIMARY_CALIBRATORS = ['URANUS','NEPTUNE','MARS'])
   if found, use first one. This sets containsPrimaryCalibrator to True or False. Do we need this? Currently we throw out tracks in main flux() program if they don't contain a primary calibrator.
7. For each antenna in refant list, try:
   1. passband calibrate (calibratePassband, passing refant array to antennasOnline variable) - only proceed if successful. Writes pb.mir
   2. measure fluxes from primary calibrator:
      1. Uvcat just primary calibrator and phase selfcal. Writes primary.mir
      2. loop over calibrators:
         1. uvcat just the calibrator and phase selfcal. Write source.mir
         2. run bootflux in each wideband window, keeping lists of freq, flux, errors, numbers of time pts:

            bootflux vis=primary.mir,source.mir primary= badres=BADRES line=chan,1,startchan,nchan taver=TAVER

            where BADRES is 30, TAVER is 10.0

         3. checks for bootflux variation of more than 25% - if so, writes just that calibrator to the fluxcal.log file with failureMode bootfluxVar
         4. check that there are a good number of time pts in each wb win
         5. print fluxes to the screen
         6. calculate spectral index and print result:
            1. if only 1 window, spec indx = -1
            2. if 2 windows, spec indx = slope, spec indx error = -1
            3. if > 2 windows, getSpecIndx and spec indx error is the fit error (added in quadrature with the rms time variation if there is more than 1 time point)
   3. if no primary calibrator, big else loop: for each calibrator create source.mir and phaseSelfcal:
      run 'fake' bootflux command to cal info - bootflux vis=source.mir primary=calname line=1,1,nchan taver=TAVER
      use bootflux output to set utsource and elsource
      set cal flux using amp from uvamp() function. I've been ignoring this part - what to do here?
   4. if success is True at the end of this, break - otherwise, continue looping, trying the next refant
      success is set to False before loop through calibrators so success=True at end of loop means that at least one calibrator flux was derived successfully. I guess this is sufficient?
8. summarize results to screen and add to database if db=True and add to our log files

def addDatabase(results)
add results to database

def checkFloat(x, ndec)
returns x converted to string with ndec numbers after the decimal place

def checkString(s)
format string for writing to database it seems

def containsFatalError(logfile)
returns True or False, searching logfile for phrase 'Fatal Error'

def createMiriadFile(vis, out, purpose=None, minfreq=None, maxfreq=None, windows=None, sources=None)
writes miriad file 'out' issuing select commands based on inputs - default is just select='-source(noise)'

def executeCommand(command, log=True, printlog=False)
runs command, returning output if log=True and printing to screen if printlog=True. Defaults are log=True and printlog=False.

def formatMiriadDate(timestamp)
converts miriad date format to postgreql

def getCorrelatorConfigurations(vis)
uses uvindex output to read frequency configurations, returning a list of configurations:

• configurations = [config1, config2, …]
• where each config is a dict of dicts:
  • config[window] is a dict containing 'freq' 'nchan' 'bw'
  • where bw is calculated from the freq. interval: bw = abs(interval * nchan) * 1000.0 (in MHz)

def getSourceList(vis)

• returns sourcelist (list of dicts of srcname and intent) and list of acceptable refants
• uses listobs to accomplish this
• determination of refants is COMPLICATED!!
  • antennasOnline = antenna list from top of Chronology of Observations table in listobs ouput
  • then read antenna locations from listobs output (E,N,U)
  • for each antenna, computes median distance to each of other antennas using locations (dist)
  • refant list is just all the online antennas sorted by median distance to other ants, in ascending order - why???

def getVarplt(mirfile, var)
returns median of the array of variable var values (by parsing the output of varplt)

  varplt vis=mirfile log=var.log yaxis=var

def initializeResults(refant=None)
initializes results dict (for bootflux output) with following parameters (set to 'null' or False) - descriptions below: (source = the calibrator that you want the flux for; cal = primary flux calibrator (ie. planets))

• results['cal'] = primary calibrator name (ie. planet)
• results['calflux'] = seems that this is not currently set
• results['elCal'] = median elevation of calibrator
• results['elSource'] = median elevation of source for 'good' fluxes
• results['flux'] = median of 'good' fluxes (see parseBootflux)
• results['freq'] = mean freq of all wideband channels (set at beginning of procesFile)
• results['method'] = planet if prim cal available, online if not (assumes gains=1)
• results['obsblock'] = obsblock name
• results['ok'] = True if fluxes are found by parseBootflux function
• results['origin'] = ie. CARMA8 for refant set to 8
• results['scale'] = median of Jy/K values from calibrator measurements
• results['rmsFlux'] = rms of 'good' fluxes (see parseBootflux)
• results['rmsScale'] = rms of Jy/K values from calibrator measurements (set to 'null' if == 0.0)
• results['skyrms'] = median rmspath value during 'source' observations (using getVarplt on source.mir)
• results['source'] = source name (calibrator you want to know the flux of)
• results['tau230'] = median tau230 value during 'source' observations (using getVarplt on source.mir)
• results['tsys'] = median of tsys values for source
• results['utCal'] = middle date/time of all cal observations
• results['utSource'] = middle date/time of source observations

def parseBootflux(bootfluxLog, refant)
parse bootflux output in bootfluxlog, filling results dict initialized with initializeResults(refant=refant): 'good' fluxes are those that have errors between 0 and 5*median(error) and have elevation ⇐ ELEVATION_MAX (80). Sets following parameters of results: cal, elCal, elSource, scale, rmsFlux, rmsScale, source, tsys, utCal, utSource

def phaseSelfcal(mirfile, refant, interval, amp=False)
does selfcal vis=mirfile refant=refant interval=interval (options=amp if amp=True)

def selectSources(sourceList, purpose, planets=True, primaryOnly=False)
Find sources that have specify intents

• sourceList : Output from getSourceList()
• purpose : Miriad purpose flag.
• planets : Include planets
• primaryOnly : Select primary calibrators only

returns sources : List of sources

def uvamp(mirfile, sourcename)
runs uvamp and reads output, returning amplitude and sigma (using just 1 bin in uvdist, so average of entire dataset)

def ampFlagPlanet(mirfile, srcname, winpar)
Peforms amplitude flagging (using calcAmpFlagging) on the planet in dataset mirfile for the window specified by winpar. Does the following steps:

1. use uvmodel options=divide to divide uv data by planet model → wb_planet.mir
2. calculate the acceptable uvdistance values: 0 to nullfrac*first null position
3. uvaver with large time interval to get 1 point per baseline → wb_planet_avg.mir
4. calcAmpFlagging allowing minimum, maximum values of 0.2, 5.0

def calcAmpFlagging(mirfile, srcname, winpar, uvindex_log=None, mirfileToFlag=None, minmax=None, winparToFlag=None, applyAbsCut=True, debug=False, dryrun=False)
General task for performing amplitude flagging for srcname, winpar in mirfile. Does the following:

1. read average amplitudes using getAvgAmps
2. check that there are > 0 non-blank values
3. define allowable values: (can be over-ridden by minmax)
   • min = median / 4.0
   • max = min(absolute cut, median * 4.0) where absolute cut = 100 Jy for 3C273, 3C279, 3C454.3, 3C84 or 50 Jy otherwise
4. for each baseline/time point, determine good or bad (within min/max)
5. Look for trends:
   • if > 50% of time slices are bad for a given baseline - mark baseline completely bad
   • if > 50% of baselines in a time slice are bad - mark entire time bad
   • if > 50% of baselines for an antenna in a time slice are bad - mark entire antenna bad in that time slice
   • if an antenna is bad in > 50% of time slices - mark entire antenna bad
6. Write flagging commands for the trends found (described above) as well as any individual baseline/time points that were marked bad.
7. Print a summary of time intervals and antennas found to be bad.
8. Return somethingleft, median flux where somethingleft is a boolean indicating whether there is enough left of this source/window to proceed. The criteria for somethingleft to be True are:
   • at least 1 time chunk remains
   • > 50% of antennas remain

def checkForPrim()
reads listobs output to find sources in the dataset. Checks both for the number of planets in the observation (given by PLANETS) and if there is a primary calibrator (PRIMARY_CALIBRATORS). Returns True if there is a primary calibrator and 0 or 1 planets in the observation. (We exclude observations with multiple planets here since the code is not equipped to handle this.)

def getAntsOnline()
returns a list of strings of the antenna numbers online for the track, derived from the listobs output - must exist!

def getAvgAmps ( vis_file, source, chanspec, uvindex_log=None, plots=False)
Using uvindex and smauvplt, finds the average amplitude on each baseline in each time chunk that source is observed in vis_file (in the channels specified by chanspec).

def getFidFreq(lofreq)
returns fiducial frequency (30,112,230 GHz for 1cm,3mm,1mm) for grade calculation

def getGrade(vis)
returns weather grade of the track (calculated a la quality script as of 9/14/11)

def getSpecIndx(nu, f, ferr)
Calculates the spectral index using chi-squared fitting of the frequencies (nu), fluxes (f) and errors in the fluxes (ferr). Assumes f = C * nu^a where C is a constant and a is the spectral index → log(f) = log(C) + a * log(nu) = b + a * log(nu). Returns list: [spec indx, spec indx err, constant, constant err]

def HMS_to_decimal(s)
returns decimal hours of date string s ('HH:MM:SS.S')

def loFreq ()
returns the LO frequency (GHz) from listobs output file (listobs.log - this must exist already!)

def maxBaseline ()
returns max basline length (m), calculated from listobs output file (listobs.log - this must exist already!). This is used in calculating grade (based on quality script)

def sh(s)
Chat's command-running thing

def sh_pipe(s)
Another of Chat's command-running things

def trackLength ()
returns track length (hours) calculated from listobs output file (listobs.log - this must exist already!)