Amber's notes on fluxes_temp.py version of 11/22/11. Please update/edit as we work on the code!

Useful information for navigating all these scripts:

• fluxes.py, fluxes_temp.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

Major changes to do:

• the various ways miriad data are written should be documented somewhere - DONE
  • AB added this to CARMA wiki here
• make script handle all dataset types properly (some of this is done with obsmode, set in flux() according to onlineAntennas - CARMA15, CARMA23 or sci2):
  • CARMA15 - standard - includes single pol, dual pol, full stokes and pacs
    • pacs - run separate, standard reductions on both files
      • sci1 track is called, ie. c0835.1A_90NGC423.1.mir
      • sci2 track is called, ie. c0835.1A_30NGC423.PACS.1.mir
      • grading - since we don't utilize pacs in reduction, I (AB) don't think we need to use special grading criteria
    • full stokes - AB added select=pol(RR,LL) to bootflux call to make sure we don't process RL and LR in full stokes datasets
  • CARMA23 - need to flag LSB (done in processFile)
    • MAXSENS - what is this? do we need to worry about it?
  • sci2 - flag LSB/USB depending on frequency (done in processFile)
• culling data that we would not want to pass to FluxSource.cat:
  • check for primary calibrator DONE 11/22 AB
  • check for grade > C, track > 2 hours (already there)
  • check for actual science track (so as not to run on test scripts) - look for filename starting with cx, c0, c1, c2 DONE 11/22 AB
• fix issue with multiple correlator configurations - DONE 11/22 AB
  • had to uvcat just wideband windows first to get around this
• error reporting
• an algorithm to interpret bootflux output?
• make code more robust:
  • flag for amplitude so that code isn't broken by 1 bad baseline - DONE 11/22 AB
  • other things??
• MARS in A array is mostly resolved - bootflux fatal error → make program catch this!
• things to change/add from Mar 2012 discussion:
  • catch bootflux crashes due to resolution issues with planets
  • full spectral index implementation
  • fancy amplitude flagging
  • need to only analyze pol(LL) (need this for parsing smauvplt output for amplitude flagging

Other stuff:

• ask Peter Teuben about uvflag bug for flagging edge channels - simplify this in flagMiriadData if possible
• calflux is not currently set by parseBootFlux?
• success evaluated well in processFile?
• raising exceptions for exit reasons - should use sys.exit('some msg') or raise SystemExit, 'some msg'
• Things to ask telecon people:
  • list of gain calibrators to accept? - allow all?

1. Just wide-band windows, pol(LL) → wb.mir
2. Remove bad windows for sci2/CARMA23 obs → wb_sci2.mir
3. Amplitude flagging creates intermediate files described below, but flags either wb.mir or wb_sci2.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
4. Passband calibration → pb.mir
5. Phase selfcal planet only → primary.mir
6. Phase selfcal each calibrator → source.mir

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. get correlator configuration of original vis file with getCorrelatorConfigurations
6. get nchan and min/max freq of correlator config, make list of wideband windows
   CRITERIA: there must be at least 1 wide-band window
7. create file of just wideband data (from calibrators.mir) selecting windows and min,max freq and only LL polarization. Writes wb.mir
8. read new correlator configurations with getCorrelatorConfigurations
9. write visfile name to temporary recording file (chat/fluxcal/fluxsource_temp_sci1/2.cat)
10. 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, 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 appears to be a list of dicts describing correlator windows (ie. 'nchan')

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(mirfile, correlator, obsmode, obsblock=None, db=False, plotResult=False, ampFlag=2)

1. 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 wb_sci2.mir
2. get updated correlator configuration: config_wb
3. get source list and calibrators list (including planets)
4. flagMiriadFile(wb.mir or wb_sci2.mir, config_wb, calibrators list)
5. extract only calibrator data → calibrators.mir
6. get updated correlator configuration (config_wb) and calculate mean freq of windows (freq)
7. get potential passband calibrator list (selecSources planets=False)
8. check for primary calibrators (PRIMARY_CALIBRATORS = ['URANUS','NEPTUNE','MARS','JUPITER','VENUS'])
   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.
9. 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 3.01

         3. check that there are a good number of time pts in each wb win
         4. print fluxes to the screen
         5. 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?
10. 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 / 10.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!)