Date - 2000.09.07 Tester - S.T. Myers (AOC) Platform - Linux (kernow) Version - Stable (1.4 build 241) Note: bugs are denoted by lines with prefix: >>>BUG and queries are denoted by lines with prefix: >>>QUERY while comments are prefixed by: >>>COMMENT: ------------------------------------------------------------------------------- GOAL: continue functionalizing of C-band polarization calibration Preliminaries: Note - to reconcile the use of TQL strings msselect for calibrater selection while imager uses both fieldid/spwid and msselect use the horrible TQL strings globally - will this work? No - setdata can AND the other selections with the msselect but setimage needs fieldid and spwid. Can these be reconciled, eg. only use TQL selection in all the set-isms? >>>COMMENT: It would be desirable to reconcile the various selections by parameters and TQL strings - maybe TQL strings are the way to go. Functions so far: # # readfits.g - Function to read in uvfits to measurement set # include 'synthesis.g' function readfits(uvfile, prefix) { dowait:=T ms_name := spaste(prefix,".ms") m:=fitstoms(msfile=ms_name,fitsfile=uvfile,readonly=T,lock=F); m.close() } # end of function readfits.g # # makedirty.g - Function to make dirty map (stokes I) # # Input measurement set = prefix.ms (eg. thems.ms) # Output files named by prefix+qual (eg. cband.run1.*) # Input cell_siz and image_siz required. # This will create and destroy its own imager object. # Returns glish record of map stats # # Version STM 2000.09.06 # function makedirty(prefix='thems',qual='run1',cell_siz,imag_siz,fldid=1,spid=1) { # Use imager to make large dirty image with single ms_name := spaste(prefix,".ms") imgr := imager(ms_name) imgr.setimage(cellx=cell_siz, celly=cell_siz, nx=imag_siz, ny=imag_siz, stokes='I', msselect=selectr) imgr.setdata(spwid=1, mode="none") imgr.weight('robust') imgr.uvrange() pre_name := spaste(prefix,qual) psf_name := spaste(pre_name,".psf") bmaj:=F; bmin:=F; bpa:=F; imgr.makeimage(type='psf',image=psf_name) imgr.fitpsf(psf=psf_name, bmaj=bmaj, bmin=bmin, bpa=bpa) img_name := spaste(pre_name,".dirty") imgr.makeimage(type='observed',image=img_name) imgr.close() imgr.done() # Get stats of dirty image img_stats := imagestats( img_name ) peak_flx := img_stats.flux peak_pix := img_stats.maxpos peak_dir := img_stats.direction img_rms := img_stats.rms print 'Peak in dirty image ',peak_flx,' Jy at ',dm.time(peak_dir.m0),dm.angle(peak_dir.m1),' at ',peak_pix print 'Map rms ',img_rms # return glish record of map stats return img_stats } # End of function makedirty.g # # imagestats.g - Function to get image stats # # Input image name required # Returns glish records with max, maxpos, maxdirec, rms # function imagestats( img_name ) { # Use image constructor to access dirty image and find peak im := image(img_name) # Search inner quarter regoff := im.shape()/2; regoff[3:4] := [1,1]; reglen := im.shape()/4; reglen[3:4] := [0,0]; reg_blc := regoff - reglen reg_trc := regoff + reglen; r := drm.box(blc=reg_blc,trc=reg_trc) im.statistics(statsout=img_stat,list=T,region=r) img_max := img_stat.maxpos img_dir := im.coordmeasures(pixel=img_max).direction img_flx := img_stat.max # Get rms noise level from an outer region regoff := im.shape()/4; regoff[3:4] := [1,1]; reglen := im.shape()/6; reglen[3:4] := [0,0]; reg_blc := regoff - reglen reg_trc := regoff + reglen; r := drm.box(blc=reg_blc,trc=reg_trc) im.statistics(statsout=noise_stat,list=T,region=r) img_rms := noise_stat.rms im.close() im.done() # Return glish record, with fluxes as quanta return [ flux=dq.quantity(img_flx,'Jy'), maxpos=img_max, direction=img_dir, rms=dq.quantity(img_rms,'Jy') ] } # End of function imagestats 1. Load up functions and run them include 'readfits.g' include 'makedirty.g' include 'imagestats.g' uvfile := 'CALS00FEB.CEDT' prefix := 'cband' qual := '.run1' readfits(uvfile,prefix) fldid := 1 spid := 1 cell_siz := '0.25arcsec' imag_siz := 256 makedirty(prefix, qual, fldid, spid, cell_siz, imag_siz) OUTPUT: > Peak in dirty image [value=0.0346042924, unit=Jy] Jy at > 09:27:03.035 +039.02.20.852 at [128 129 1 1] > Map rms [value=0.00541110057, unit=Jy] > [flux=[value=0.0346042924, unit=Jy], maxpos=[128 129 1 1] , > direction=[type=direction, refer=J2000, m1=[value=0.681361278, unit=rad], > m0=[unit=rad, value=2.47422495]], rms=[value=0.00541110057, unit=Jy]] 2. Lets try to get selfcal to work starting with point model I want to be able to convert the fieldid and spwid selection parameters into TQL strings, eg. selectr := sprintf('FIELD_ID==%d SPW_ID==%d',fldid,spid) print selectr > FIELD_ID==1 SPW_ID==1 I need a list of the keywords, I know FIELD_ID is one, but what is the spwid equivalent? There should be a table of these easily found in the documentation. >>>COMMENT: Should be a list of TQL keywords accessible from the documentation entries for imagr.setdata, calibrater.setdata, etc. I will use the table.summary function instead tb := table("cband.ms") tb.summary() OUTPUT: Table keywords: [ANTENNA=Table: cband.ms/ANTENNA, ARRAY=Table: cband.ms/ARRAY, FEED=Table: cband.ms/FEED, FIELD=Table: cband.ms/FIELD, OBSERVATION=Table: cband.ms/OBSERVATION, OBS_LOG=Table: cband.ms/OBS_LOG, SOURCE=Table: cband.ms/SOURCE, SPECTRAL_WINDOW=Table: cband.ms/SPECTRAL_WINDOW, SYSCAL=Table: cband.ms/SYSCAL, WEATHER=Table: cband.ms/WEATHER, SORTED_TABLE=Table: cband.ms/SORTED_TABLE, SORT_COLUMNS=FIELD_ID ARRAY_ID SPECTRAL_WINDOW_ID TIME] Columns: UVW ANTENNA1 ANTENNA2 ARRAY_ID CORRELATOR_ID EXPOSURE FEED1 FEED2 FIELD_ID FLAG_ROW INTERVAL OBSERVATION_ID PULSAR_BIN PULSAR_ID SCAN_NUMBER SPECTRAL_WINDOW_ID TIME WEIGHT DATA SIGMA WEIGHT_SPECTRUM FLAG MODEL_DATA CORRECTED_DATA IMAGING_WEIGHT UVW keywords: [UNIT=m, MEASURE_TYPE=UVW, MEASURE_REFERENCE=] ANTENNA1 keywords: [UNIT=, MEASURE_TYPE=] ANTENNA2 keywords: [UNIT=, MEASURE_TYPE=] ARRAY_ID keywords: [UNIT=, MEASURE_TYPE=] CORRELATOR_ID keywords: [UNIT=, MEASURE_TYPE=] EXPOSURE keywords: [UNIT=s, MEASURE_TYPE=] FEED1 keywords: [UNIT=, MEASURE_TYPE=] FEED2 keywords: [UNIT=, MEASURE_TYPE=] FIELD_ID keywords: [UNIT=, MEASURE_TYPE=] FLAG_ROW keywords: [UNIT=, MEASURE_TYPE=] INTERVAL keywords: [UNIT=s, MEASURE_TYPE=] OBSERVATION_ID keywords: [UNIT=, MEASURE_TYPE=] PULSAR_BIN keywords: [UNIT=, MEASURE_TYPE=] PULSAR_ID keywords: [UNIT=, MEASURE_TYPE=] SCAN_NUMBER keywords: [UNIT=, MEASURE_TYPE=] SPECTRAL_WINDOW_ID keywords: [UNIT=, MEASURE_TYPE=] TIME keywords: [UNIT=s, MEASURE_TYPE=EPOCH, MEASURE_REFERENCE=TAI] WEIGHT keywords: [UNIT=, MEASURE_TYPE=] DATA keywords: [UNIT=, MEASURE_TYPE=] SIGMA keywords: [UNIT=, MEASURE_TYPE=] WEIGHT_SPECTRUM keywords: [UNIT=, MEASURE_TYPE=] FLAG keywords: [UNIT=, MEASURE_TYPE=] MODEL_DATA keywords: [CHANNEL_SELECTION=[[1:2,] 0 0 1 1]] tb.close() Ah - the keyword I need is SPECTRAL_WINDOW_ID selectr := sprintf('FIELD_ID==%d && SPECTRAL_WINDOW_ID==%d',fldid,spid) OK, lets try the following function: # # pointcal.g - Function to do selfcal on data # # Makes calibrater object using point source with known flux # Returns glish record of last residual map stats # # Note input cal_flux is vector [ I, Q, U, V ] - default=[1.0,0,0,0] # # Currently works off of field and spwid params and constructs the TQL # selection string. # include 'mathematics.g' include 'gainstats.g' function pointcal(prefix='thems',qual='run1',fldid=1,spid=1,cal_flux=[1.0,0,0,0]) { # Some local parameters ms_name := spaste(prefix,".ms") pre_name := spaste(prefix,qual) sn_name := spaste(pre_name,".sn") # Data selection parameters selectr := sprintf('FIELD_ID==%d && SPECTRAL_WINDOW_ID==%d',fldid,spid) # Make imager object imgr:=imager(ms_name) # put the model visibilities in the MODEL_DATA column print 'Using point source of flux ',cal_flux imgr.setjy(fieldid=fldid,spwid=spid,fluxdensity=cal_flux) imgr.close() # Make calibrater object cl:=calibrater(ms_name) cl.setdata(mode="none", nchan=1, start=1, step=1, mstart=[value=0.0, unit="km/s" ], mstep=[value=0.0, unit="km/s" ], uvrange=0, msselect=selectr); # Set selfcal parameters for gain solution, then solve cl.setsolve (type="G", t=60.0, phaseonly=F, table=sn_name, append=F); cl.solve() print 'Made calibration table ',sn_name # Apply solution to data cl.setapply (type="G", t=60.0, table=sn_name, msselect=selectr); cl.correct() # Done with calibrater cl.close() # Get the statistics g_stat := gainstats(sn_name) # return glish record of gain cal stats return g_stat } # end of function pointcal.g >>>COMMENT: This took a long time to debug due to problems with the selection string - I forgot to put the Boolean && in between the selection elements. Ger's TaQL description was confusing, though it did mention Boolean operators. The example under the URG entry for imagr.setdata did have an example with two selections && together. The examples in calibrater.setdata and .setapply should do the same. The stats are computed in the function # # gainstats.g - Function to report back calibration gain statistics # include 'mathematics.g' function gainstats(sn_name) { # open the table gtable := table(sn_name) print 'Opened calibration table ',sn_name # now check through the calibration table to extract the good gains gains := gtable.getcolslice(columnname='GAIN',blc=[1,1,1],trc=[1,1,2]) goods := gtable.getcol(columnname='SOLUTION_OK') # done with table gtable.done() # mask out good gain sols into vector, do stats ggain := gains[goods] ngain := len(ggain) amps := real( sqrt( ggain*conj(ggain) ) ) amprms := stddev( amps ) ampmin := min( amps ) ampmax := max( amps ) phases := arg( ggain )*180.0/pi pharms := stddev( phases ) phamin := min( phases ) phamax := max( phases ) # return glish record of phase cal stats return [ ampstats=[ rms=amprms, min=ampmin, max=ampmax], phastats=[rms=dq.quantity(pharms,'deg'), min=dq.quantity(phamin,'deg'), max=dq.quantity(phamax,'deg')], number=ngain ] } # end of function gainstats.g We run these with include 'pointcal.g' # Set flux of 0927+390 (fieldid=1) to 10.65 Jy as determined using AIPS # and posted at http://www.aoc.nrao.edu/~smyers/calibration/C_band.html # Set Q and U (and V) fluxes to 0.0 initially cal_flux := [10.65,0.0,0.0,0.0] pointcal(prefix, qual, fldid, spid, cal_flux) Message Log for solutions: Starting calibrater::solve Initializing solvable electronic gain terms (G-matrix) For interval of 60 seconds, found 2 slots Solving for G G Jones Slot=1, 0927+390, spw=1: 24-Feb-2000/09:43:00 to 24-Feb-2000/09:43:59 G Jones Initial fit per unit weight = 1.48921 Jy, sum of weights = 15500 G Jones Final fit per unit weight = 0.0225948 Jy after 9 iterations G Jones Slot=2, 0927+390, spw=1: 24-Feb-2000/09:44:00 to 24-Feb-2000/09:44:59 G Jones Initial fit per unit weight = 1.50231 Jy, sum of weights = 3700 G Jones Final fit per unit weight = 0.0230248 Jy after 9 iterations Storing G matrix in table cband.run1.sn Finished calibrater::solve Function output: (edited for clarity) > [ampstats=[rms=0.332134518, min=0.292028814, max=1], > phastats=[rms=[value=41.1757172, unit=deg], > min=[value=-136.266316, unit=deg], > max=[value=127.446219, unit=deg]], > number=108] 3. Things that worry me: the max gain amplitude of 1, the two slot outputs from solve() - there should have been more if these are time intervals. Things that dont worry me: the high phase rms (41deg) on 60s intervals! Note that we have a full Jones matrix solution, should grab those from the .sn table when I eventually do the poln solutions. Use GUI table browser to look at cband.run1.sn, examine one of the gain entries, eg. Antenna 12 : Gain[1,1] = -0.332683682i Gain[2,2] = -0.0166011564+0.337170303i Gain[1,1],Gain[2,1] = 0+0i However, AN1 (row 2) has the gains: Antenna 2: Gain[1,1] = 12+0.1206434i Gain[2,2] = -0.236283958+0.240442574i How does the 12 jibe with the reported stats? Note all the other entries have purely imaginary Gain[1,1]. Funny. >>>BUG: when I invoked the .close method and then did Done on the GUI, got the following droppings on the glish window: "newab.g", line 500: warning, operand to .cleanup is not a record "newab.g", line 500: error, F is not a function value NEXT: Lets find those missing time slots! And try visplot to look at results. Is there an equivalent of SNPLT to graphically view the solutions?