Date: Wed, 7 May 2008 10:42:38 -0600 From: Debra Shepherd Subject: Re: CASA mm tutorial This was a BIMA SONG galaxy. Michele Thornley gave it to me in miriad format for the 2nd offline user test. Then Joe made it into a regression script after the user test. >From my notes: N4826 - BIMA SONG Data 16apr98 source=ngc4826 phasecal=1310+323 fluxcal=3c273, Flux = 23 Jy on 16apr98 passcal= none - apparently, this is considered optional... odd 22apr98 source=ngc4826 phasecal=1230+123 fluxcal=3c273, Flux = 23 Jy on 22apr98 passcal= none # miriad: source Vlsr = 408; delta is 20 km/s ? NOTE: This data has been filled into MIRIAD, line-length correction done, and then exported as separate files for each source. 3c273 was not line length corrected since it was observed for such a short amount of time that it did not need it. The date were filled into MS format with the mirfiller tool in aips++ and then concatenated. I have the original data. Some editing was needed. Again from my original notes: plot := msplot(msfile='all.ms', edit=T); # field 1, spwid 1 (one wide-band channel) # no obvious bad data, BL3-7 is low but not excessive. # field 2, spwid 2 (one wide-band channel) # again, nothing obviously bad # fields 3,4,5,6,7,8,9, spwids 3,4,5,6 # this data has NOT been pre-edited. - GOOD. # - monster bad correlator glitch: BL 3-9, field 9 # timerange: 1998/04/16/06:19:00 to 1998/04/16/06:21:00 # - need to flag first and last 2 channels of source data # (no band pass calibration) I flagged the correlator glitch and the first and last channels and my notes say the data look 'nicely behaved.' ########################################################### ## Mosaic field spacing looks like: ## ## F3 (fieldid 4) F2 (fieldid 3) ## ## F4 (fieldid 5) F0 (fieldid 1) F1 (fieldid 2) ## ## F5 (fieldid 6) F6 (fieldid 7) ## ## 4x64 channels = 256 channels (any overlap?) ## Image central 192 channels, 4 chan averaging => 48 chans ## Start on chan 32 ## ## Primary Beam should be about 1.6' FWHM (7m dishes, 2.7mm wavelength) ## Resolution should be about 5-8" -- Date: Wed, 7 May 2008 16:10:56 -0600 From: Debra Shepherd Subject: Re: CASA mm tutorial Hi Steve, I have taken the miriad line-length corrected data for one day of the NGC 4826 observations, used miriad to write it to fits files, then used casa to import them. I looked at the data in plotxy and identified what is bad - it looks like a good dataset to practice on. You can easily use this as a guide to process the first day (apr16) if this one turns out to not be suitable. Below I give the script in miriad I used and the casa commands I used to fill them and then concatenate all the ms's into a single ms that should be good to go for the regression script you have. I came across 3 problems though - 1. I have compared the frequencies at every stage of the convert and fill process and when I fill the data into CASA I get a several MHz frequency shift from what my notes report was obtained in AIPS++. I cc Brian because of this (we found a similar frequency shift in ATF spectral line observations of Orion that Al Wootten took a few weeks back). Please look at my script and all the frequency reporting I did - I may have made a mistake. If you concur however, we have a problem with the frequency determination in CASA fill. 2. Plotants task gives a severe error even though I get a plot saying the name of the telescope is not in the MS - but it was there in the fits file so this looks like a fits filler bug. I don't know what this will do for subsequent processing. 3. The plotants task chokes - it plots an extra point in never-never land (far from the array central). Attached is my detailed script (located at: /home/sola/dss/support/summer.school.08/bima.song.ngc4826/map The data I used for this test are in the directory: /home/sola/dss/support/summer.school.08/bima.song.ngc4826/data.98apr22 -------------------------------------- Cheers, Debra ======================================================================= Debra's Script and Notes ======================================================================= 7may08 - NRAO Synthesis Summer School CASA tutorial for millimeter reduction and mosaic imaging N4826 - BIMA SONG Data (courtesy of Michele Thornley) 16apr98 source=ngc4826 phasecal=1310+323 fluxcal=3c273, Flux = 23 Jy on 16apr98 passcal= none - apparently, this is considered optional... odd 22apr98 source=ngc4826 phasecal=1230+123 fluxcal=3c273, Flux = 23 Jy on 22apr98 passcal= none source Vlsr = 408; line width is 20 km/s NOTE: This data has been filled into MIRIAD, line-length correction done, and then exported as separate files for each source. 3c273 was not line length corrected since it was observed for such a short amount of time that it did not need it. /home/sola/dss/support/summer.school.08/bima.song.ngc4826/data.98apr22 # miriad files are located in 2 directories (one for each day): 98apr16/ has 1310+323.ll/ 3c273/ ngc4826.ll/ 98apr22/ has 1230+123.ll/ 3c273/ ngc4826.ll/ Just work with the 98apr22 data (arbritrary choice): --------------------------------------- uvlist vis=3c273 options=full,var nschan : 32 32 32 32 32 32 32 32 nschan : 32 32 32 32 32 32 32 32 sfreq : 112.07 111.97 111.87 111.77 111.67 sfreq : 111.57 111.47 111.37 114.82 114.92 sfreq : 115.02 115.12 115.22 115.32 115.42 sfreq : 115.52 uvlist vis=1230+123.ll options=full,var nschan : 32 32 32 32 32 32 32 32 nschan : 32 32 32 32 32 32 32 32 sfreq : 112.07 111.97 111.87 111.77 111.67 sfreq : 111.57 111.47 111.37 114.82 114.92 sfreq : 115.02 115.12 115.22 115.32 115.42 sfreq : 115.52 uvlist vis=ngc4826.ll options=full,var nschan : 64 64 64 64 64 64 64 64 sfreq : 111.94 111.85 111.76 111.67 114.95 sfreq : 115.04 115.13 115.22 # For the calibrators: I need a single wide band data in the USB (where CO line is in source) There are 8 subbands in the USB, each with 32 channels - average these all together to form a single wide band continuum channel that can be used for gain and flux calibration. #### NOTE: spectral windows of calibrators appear to be identical. # For the source: I need the CO line in the USB, It should be broad, covering all four spectral windows (5,6,7,8) # in miriad, convert these miriad format files into fits files: fits in=3c273 op=uvout out=3c273.fits line=channel,1,1,256 select='window(9,10,11,12,13,14,15,16)' Polarisations copied: YY. 2025 visibilities copied to 2025 output records. less 3c273.fits gives: CTYPE3 = 'STOKES ' / CRVAL4 = 1.15213996254E+11 / CRPIX4 = 1.00000000000E+00 / CDELT4 = 8.00000000000E+08 / fits in=1230+123.ll op=uvout out=1230+123.fits line=channel,1,1,256 select='window(9,10,11,12,13,14,15,16)' Polarisations copied: YY. 15480 visibilities copied to 15480 output records. less 1230+123.fits gives: CTYPE3 = 'STOKES ' / CRVAL4 = 1.15214060656E+11 / - this is now slightly different, must use frequency tolerance CRPIX4 = 1.00000000000E+00 / CDELT4 = 8.00000000000E+08 / fits in=ngc4826.ll op=uvout out=ngc4826-1.fits line=channel,64,1 select='window(5)' Polarisations copied: YY. 13455 visibilities copied to 13455 output records. fits in=ngc4826.ll op=uvout out=ngc4826-2.fits line=channel,64,1 select='window(6)' fits in=ngc4826.ll op=uvout out=ngc4826-3.fits line=channel,64,1 select='window(7)' fits in=ngc4826.ll op=uvout out=ngc4826-4.fits line=channel,64,1 select='window(8)' casapy # Start CASA software ########################################################### ## Fill Data just from 22apr98: ## importuvfits(fitsfile='3c273.fits', vis='3c273.ms') listobs(vis='3c273.ms') Fields: 1 ID Code Name Right Ascension Declination Epoch 0 3C273 12:29:06.70 +02.03.08.60 J2000 SpwID #Chans Frame Ch1(MHz) ChanWid(kHz)TotBW(kHz) Ref(MHz) Corrs Spectral Windows: (1 unique spectral windows, 1 unique polarization setups) 0 1 LSRK 115213.996 800000 800000.001 115213.996 YY importuvfits(fitsfile='1230+123.fits', vis='1230+123.ms') listobs(vis='1230+123.ms') Fields: 1 ID Code Name Right Ascension Declination Epoch 0 1230+123 12:30:49.42 +12.23.28.04 J2000 Spectral Windows: (1 unique spectral windows, 1 unique polarization setups) SpwID #Chans Frame Ch1(MHz) ChanWid(kHz)TotBW(kHz) Ref(MHz) Corrs 0 1 LSRK 115214.061 800000 800000.001 115214.061 YY importuvfits(fitsfile='ngc4826-1.fits', vis='ngc4826-1.ms') importuvfits(fitsfile='ngc4826-2.fits', vis='ngc4826-2.ms') importuvfits(fitsfile='ngc4826-3.fits', vis='ngc4826-3.ms') importuvfits(fitsfile='ngc4826-4.fits', vis='ngc4826-4.ms') listobs(vis='ngc4826-1.ms') listobs(vis='ngc4826-2.ms') listobs(vis='ngc4826-3.ms') listobs(vis='ngc4826-4.ms') ID #Chans Frame Ch1(MHz) Resoln(kHz) TotBW(kHz) Ref(MHz) Corrs 0 64 LSRK 114948.307 1562.5 1562.564 114948.307 YY 0 64 LSRK 115038.321 1562.5 1562.564 115038.321 YY 0 64 LSRK 115128.059 1562.5 1562.564 115128.059 YY 0 64 LSRK 115218.073 1562.5 1562.564 115218.073 YY #### humm, all several MHz off from what I got before with AIPS++: ID #Chans Frame Ch1(MHz) Resoln(kHz) TotBW(kHz) Ref(MHz) Corrs 1 1 LSRD 115219.429 800000.072 800000.072 115271.2 YY 2 64 LSRD 114950.382 1562.5 100000 115271.2 YY 3 64 LSRD 115040.397 1562.5 100000 115271.2 YY 4 64 LSRD 115130.137 1562.5 100000 115271.2 YY 5 64 LSRD 115220.153 1562.5 100000 115271.2 YY ########################################################### ## Concatenate all ms's ## !cp -r 3c273.ms all.ms concat(vis='all.ms', concatvis='1230+123.ms', freqtol='10MHz') listobs(vis='all.ms') SpwID #Chans Frame Ch1(MHz) ChanWid(kHz)TotBW(kHz) Ref(MHz) Corrs 0 1 LSRK 115213.996 800000 800000.001 115213.996 YY concat(vis='all.ms', concatvis='ngc4826-1.ms') concat(vis='all.ms', concatvis='ngc4826-2.ms') concat(vis='all.ms', concatvis='ngc4826-3.ms') concat(vis='all.ms', concatvis='ngc4826-4.ms') listobs(vis='all.ms') Fields: 9 ID Code Name Right Ascension Declination Epoch 0 3C273 12:29:06.70 +02.03.08.60 J2000 1 1230+123 12:30:49.42 +12.23.28.04 J2000 2 NGC4826 12:56:44.24 +21.41.05.10 J2000 3 NGC4826 12:56:41.08 +21.41.05.10 J2000 4 NGC4826 12:56:42.66 +21.41.43.20 J2000 5 NGC4826 12:56:45.82 +21.41.43.20 J2000 6 NGC4826 12:56:47.39 +21.41.05.10 J2000 7 NGC4826 12:56:45.82 +21.40.27.00 J2000 8 NGC4826 12:56:42.66 +21.40.27.00 J2000 Spectral Windows: (5 unique spectral windows and 1 unique polarization setups) SpwID #Chans Frame Ch1(MHz) ChanWid(kHz)TotBW(kHz) Ref(MHz) Corrs 0 1 LSRK 115213.996 800000 800000.001 115213.996 YY 1 64 LSRK 114948.307 1562.5 1562.564 114948.307 YY 2 64 LSRK 115038.321 1562.5 1562.564 115038.321 YY 3 64 LSRK 115128.059 1562.5 1562.564 115128.059 YY 4 64 LSRK 115218.073 1562.5 1562.564 115218.073 YY What I had before with AIPS++: Fields: 9 ID Name Right Ascension Declination Epoch 1 3C273-F0 12:29:06.70 +02.03.08.60 J2000 fcal 2 1230+123-F0 12:30:49.42 +12.23.28.04 J2000 gcal 3 NGC4826-F0 12:56:44.24 +21.41.05.10 J2000 src 4 NGC4826-F1 12:56:41.08 +21.41.05.10 J2000 " 5 NGC4826-F2 12:56:42.66 +21.41.43.20 J2000 " 6 NGC4826-F3 12:56:45.82 +21.41.43.20 J2000 " 7 NGC4826-F4 12:56:47.40 +21.41.05.10 J2000 " 8 NGC4826-F5 12:56:45.82 +21.40.27.00 J2000 " 9 NGC4826-F6 12:56:42.66 +21.40.27.00 J2000 " Data descriptions: 5 (5 spectral windows and 1 polarization setups) ID #Chans Frame Ch1(MHz) Resoln(kHz) TotBW(kHz) Ref(MHz) Corrs 1 1 LSRD 115219.429 800000.072 800000.072 115271.2 YY 2 64 LSRD 114950.382 1562.5 100000 115271.2 YY 3 64 LSRD 115040.397 1562.5 100000 115271.2 YY 4 64 LSRD 115130.137 1562.5 100000 115271.2 YY 5 64 LSRD 115220.153 1562.5 100000 115271.2 YY ########################################################### ## Manually plot the data for visual inspection: ## plotxy(vis='all.ms'); # flux cal: field 0, spwid 0 (one wide-band channel) # ANT9 is low but not excessive between 5h10m25s and 5h10m41s # gain cal: field 1, spwid 0 (one wide-band channel) # 3rd set of integrations are bonkers and then it looks like it was flagged. # this will have to go along with source observations on either side. # source mosaic fields 2,3,4,5,6,7,8 spwids 1,2,3,4 # poor/bad data between 04:36 and 04:57 - # looks like ant 8 went slightly crazy (amps double than at # other timestamps and then, short pause while amps go back to near # normal but ant 8 still slightly high. # Next source block shows ant8 normal. # There are also some zero-amplitude points on the source data only # and these appear to have zero-phase (or some other points have # zero-phase). So these should be deleted also. This is a good dataset for identifying bad data. You can also see the noise increase as the source sets and you have to recognize that it is not bad, just noisy. # also, delete 1st 2 and last 2 end channels in source data plotants(vis='all.ms') SEVERE Error: Unable to plot telescope array the name of the telescope is not stored in the measurement set. INFO2 Number of points being plotted : 12 ############## BUG: Looks like the fill process dropped the antenna name. In the fits file it was listed as: TELESCOP= 'HATCREEK' I get the plot but did we loose the name of the BIMA array in casa??? ############## BUG: ### Interestingly enough - there are only 10 antennas but there is labeled as #12 way off in the corner that doesn't show up in listobs output. Looks like a bug. # Use central antenna number 4 for refant for future processing. ########################################################### ########################################################### ########################################################### STOP - the rest of this stuff is my old AIPS++ script that I developed for the second ALMA offline test - likely the regression test steps are somewhat different. ########################################################### ########################################################### ########################################################### ## Set the flux density of 3c273 (field 1) at 23 Jy flux density: ## imgr:=imager(filename='all.ms'); # size goes from 32 MB to 99 MB imgr.setjy(fieldid=1,fluxdensity=23.0); # 3C273-F0 spwid= 1 [I=23, Q=0, U=0, V=0] Jy, (user-specified) imgr.setdata(fieldid=1); imgr.plotvis(type='model'); imgr.done(); ########################################################### ## Manually plot the data for visual inspection: ## plot := msplot(msfile='all.ms', edit=T); # field 1, spwid 1 (one wide-band channel) # ANT9 is low but not excessive between 5h10.5m and 5h11m # field 2, spwid 1 (one wide-band channel) # again, nothing obviously bad - src is setting after hr 10:30 # fields 3,4,5,6,7,8,9, spwids 2,3,4,5 # poor/bad data between 04:36 and 04:57 - # looks like ant 8 went slightly crazy (amps double that at # other timestamps and then, short pause while amps go back to near # normal but ant 8 still slightly high. # Next source block shows ant8 normal. gcal normal in this range. # just delete source data # antenna plot shows this is the same configuration as the 16apr98 data # Use central antenna number 4 for refant # also, delete 1st 2 and last 2 end channels af:=autoflag(msname='all.ms'); af.setdata(); af.setselect(ant=8,timerng="23-Apr-1998/04:36:00 23-Apr-1998/04:57:00"); af.run(); af.resetall(); af.setselect(chan=[1,2],spwid=[2,3,4,5]); af.run(); af.resetall(); af.setselect(chan=[63,64],spwid=[2,3,4,5]); af.run(); af.done(); # too many logger messages... plot := msplot(msfile='all.ms', edit=f); # OK, looks reasonable. ########################################################### ########################################################### ## Derive gain calibration solutions, try VLA-like calibration: ## cal:=calibrater(filename='all.ms'); cal.reset(); cal.setdata(msselect='FIELD_ID IN [1,2] && SPECTRAL_WINDOW_ID IN [1]'); cal.setsolve(type='G', t=0, refant=4, table='all.gcal'); cal.solve(); # Found 14 good G Jones solutions. cal.plotcal(tablename='all.gcal',plottype='AMP', spwids=1); cal.plotcal(tablename='all.gcal',plottype='PHASE',spwids=1); # OK ########################################################### ## No bandpass calibration - done on line... ## ########################################################### ## Transfer the flux density scale: ## Found a new parameter in weekly RefMan documentation: refspwmap ## cal.fluxscale(tablein='all.gcal', tableout='all.fcal', reference='3C273-F0', transfer=['1230+123-F0']); # Flux density for 1230+123-F0 in SpW=1 is: # 3.27734 +/- 0.0351508 (SNR = 93.2366) cal.plotcal(tablename='all.fcal', plottype='AMP', multiplot=F); cal.plotcal(tablename='all.fcal', plottype='PHASE', multiplot=F); ########################################################### ## Correct the source and gcal data: ## Use new parm spwmap to apply gain solutions derived from spwid1 ## to all other spwids... cal.reset(); cal.setdata(msselect='FIELD_ID IN [2] && SPECTRAL_WINDOW_ID in [1]'); cal.setapply(type='G', t=0.0, table='all.fcal'); cal.correct(); cal.reset(); cal.setdata(msselect='FIELD_ID IN [3:9] && SPECTRAL_WINDOW_ID in [2:5]'); cal.setapply(type='G', t=0.0, table='all.fcal', spwmap=[1,1,1,1,1,1]); cal.correct(); # Applying G Jones calibration from spw=1 to spw=2 # Applying G Jones calibration from spw=1 to spw=3 # Applying G Jones calibration from spw=1 to spw=4 # Applying G Jones calibration from spw=1 to spw=5 cal.done(); ########################################################### ## Look at calibrated source data: ## plot := msplot(msfile='all.ms', edit=T); # source really scraping bottom at the end, clip source above 180 Jy # still have variable weather... # Summary shows that BIMA correlations are 'YY' - af:=autoflag(msname='all.ms'); af.setdata(); cliprec['ABS YY'] := [0,180] af.setselect(clip=cliprec); af.run(); af.resetall(); af.done(); # there are 322 'chunks' of data, each with a flagging report - too many # to review, can't tell what happened. I see no flags during my spot check. plot := msplot(msfile='all.ms', edit=T); # Problem opening the measurement set. The error was: # F File: note.g, Line 59 Stack: throw(), msplot.g line 6077 .(), msplot.g line 6115 msplot() plot := msplot(msfile='all.ms', edit=T); # Problem opening the measurement set. The error was: # Cannot create the flag table called[i_am_unset=i_am_unset]:create: # Table all.ms.flags.1 already open; maybe close it first # Looks like autoflag and msplot are not getting along. # exit aips++ and restart plot := msplot(msfile='all.ms', edit=T); # the data is not clipped ... # All documentation says to use a min of 1e-6 instead of 0 - maybe there is # a hidden bug that won't allow you to use 0 as the min and the # doc just gives the workaround without mentioning this problems. # Check: af:=autoflag(msname='all.ms'); af.setdata(); cliprec['ABS YY'] := [1e-6,180] af.setselect(clip=cliprec); af.run(); af.resetall(); af.done(); plot := msplot(msfile='all.ms', edit=T); # NO flagging was done. But at least msplot came up after autoflag # Maybe the YY specification is not being selected correctly: af:=autoflag(msname='all.ms'); af.setdata(); cliprec['ABS XX'] := [1e-6,180] af.setselect(clip=cliprec); *clip = [ABS YY=[1e-06 180] , ABS XX=[1e-06 180] ] # select won't take straight XX specification, no resets help # see what happens: af.run(); # LocalExec::SetStatus: abnormal child termination for # /home/ballista/aips++/stable/linux_rh9/bin/autoflag # Server 'autoflag' has failed unexpectedly! # You will need to create the relevant tool again. # If that causes unexpected behavior, please restart AIPS++ # Please submit a bug-report using bug() if you can reproduce the problem. # : Server 'autoflag' has failed unexpectedly! # You will need to create the relevant tool again. # If that causes unexpected behavior, please restart AIPS++ # Please submit a bug-report using bug() if you can reproduce the problem. # File: servers.g, Line 1021 # Stack: .(), autoflag.g line 436 # .() # stupid error message. Doesn't help me figure out what goes. # submit bug report and clip in msplot again. plot := msplot(msfile='all.ms', edit=T); # clip source above 180 Jy, clip gcal above 17Jy. ########################################################### ## Split out calibrated target source and gcal data: ## mset:=ms(filename="all.ms", readonly=F); mset.split(outputms="gcal.23apr.ms", fieldids=2, spwids=1, nchan=1, start=1, step=1, whichcol="CORRECTED_DATA"); mset.done(); mset:=ms(filename="all.ms", readonly=F); mset.split(outputms="src.23apr.ms", fieldids=[3,4,5,6,7,8,9], spwids=[2,3,4,5], nchan=64, start=1, step=1, whichcol="CORRECTED_DATA"); mset.done(); ########################################################### ## Get quick summary with default catalog tool: ## dc.summary('gcal.23apr.ms'); # Fields: 1 # ID Name Right Ascension Declination Epoch # 1 1230+123-F0 12:30:49.42 +12.23.28.04 J2000 # Data descriptions: 1 (1 spectral windows and 1 polarization setups) # ID #Chans Frame Ch1(MHz) Resoln(kHz) TotBW(kHz) Ref(MHz) Corrs # 1 1 LSRD 115219.429 800000.072 800000.072 115271.2 YY dc.summary('src.23apr.ms'); # Fields: 7 # ID Name Right Ascension Declination Epoch # 1 NGC4826-F0 12:56:44.24 +21.41.05.10 J2000 # 2 NGC4826-F1 12:56:41.08 +21.41.05.10 J2000 # 3 NGC4826-F2 12:56:42.66 +21.41.43.20 J2000 # 4 NGC4826-F3 12:56:45.82 +21.41.43.20 J2000 # 5 NGC4826-F4 12:56:47.40 +21.41.05.10 J2000 # 6 NGC4826-F5 12:56:45.82 +21.40.27.00 J2000 # 7 NGC4826-F6 12:56:42.66 +21.40.27.00 J2000 # Data descriptions: 4 (4 spectral windows and 1 polarization setups) # ID #Chans Frame Ch1(MHz) Resoln(kHz) TotBW(kHz) Ref(MHz) Corrs # 1 64 LSRD 114950.382 1562.5 100000 115271.2 YY # 2 64 LSRD 115040.397 1562.5 100000 115271.2 YY # 3 64 LSRD 115130.137 1562.5 100000 115271.2 YY # 4 64 LSRD 115220.153 1562.5 100000 115271.2 YY ########################################################### ## Image: ## ------ ## include 'imager.g' imgr:=imager(filename='gcal.23apr.ms'); imgr.setdata(fieldid=[1],spwid=[1]); imgr.setimage(nx=256, ny=256, cellx='1arcsec', celly='1arcsec', stokes='I'); imgr.clean(algorithm="clark" , niter=100, gain=0.1, threshold='0.01Jy', model="gcal.cln" , mask='gcal.mask', image="gcal.cm" , residual="gcal.resid", interactive=T, npercycle=30); # Beam used in restoration: 8.50431 by 5.65108 (arcsec) at pa -172.211 (deg) # Changed mask each iteration - looks like this gcal is resolved big time... # Glish errors at the end as usual imgr.done(); dv.gui(); # - Not a very good image - cute little bipolar jet in this source... # gain calibration not great since this is not a point source. # Peak = 3.262 Jy - OK. # RMS = 47 mJy/bm ! horrid # this is a poor image, either there is a bad antenna in here # or the gain cal was rotten because the source is resolved. # jet brightness = 0.54Jy ~ 15% level. # DEAL WITH THIS LATER. PROCEED TO MOSAIC IMAGING. ########################################################### ## Mosaic field spacing looks like: ## ## F3 (fieldid 4) F2 (fieldid 3) ## ## F4 (fieldid 5) F0 (fieldid 1) F1 (fieldid 2) ## ## F5 (fieldid 6) F6 (fieldid 7) ## ## 4x64 channels = 256 channels ## ## Primary Beam should be about 1.6' FWHM (7m dishes, 2.7mm wavelength) ## Resolution should be about 5-8" #################################################################### dell -r src.all.ms cp -r ../16apr98/src.16apr.ms/ src.all.ms mset:=ms(filename='src.all.ms',readonly=F); # fields in different days have slightly different positions and frequencies. # In order to get the fields to have the same ID, must use # freq and direction tolerances: mset.concatenate(msfile='src.23apr.ms', freqtol='1MHz', dirtol='0.01arcsec'); mset.done(); dc.summary('src.all.ms') Fields: 7 ID Name Right Ascension Declination Epoch 1 NGC4826-F0 12:56:44.24 +21.41.05.10 J2000 2 NGC4826-F1 12:56:41.08 +21.41.05.10 J2000 3 NGC4826-F2 12:56:42.66 +21.41.43.20 J2000 4 NGC4826-F3 12:56:45.82 +21.41.43.20 J2000 5 NGC4826-F4 12:56:47.40 +21.41.05.10 J2000 6 NGC4826-F5 12:56:45.82 +21.40.27.00 J2000 7 NGC4826-F6 12:56:42.66 +21.40.27.00 J2000 Data descriptions: 4 (4 spectral windows and 1 polarization setups) ID #Chans Frame Ch1(MHz) Resoln(kHz) TotBW(kHz) Ref(MHz) Corrs 1 64 LSRD 114950.382 1562.5 100000 115271.2 YY 2 64 LSRD 115040.397 1562.5 100000 115271.2 YY 3 64 LSRD 115130.137 1562.5 100000 115271.2 YY 4 64 LSRD 115220.153 1562.5 100000 115271.2 YY - spwid mapping OK now Make all calculations myself (aips++ won't help me out here). - each chan = 1.5625 MHz wide, 64 channels = 100 MHz 1st chans are separated by 90.015 MHz ==> 5.76 channel overlap Is the 5.76 because the above listing only gives freq at chan center? Assume a 6 chan overlap Approx. Total number of channels = 58+52+52+58=220 channels... Drop the first and last 6 anyway - so have 208 useful channels # Verify velocity coverage, total number of channels, and spwid overlap, # no mosaic NOTE: if you don't use channels in spwid 4, you cannot select spwid 4 in setdata/setimage or you get SEVERE errors. Kumar was supposed to put in a warning on this - see if it is there. No, same error as before. Only select spwid 1-3 imgr:=imager(filename='src.all.ms'); imgr.setdata(mode='channel', fieldid=[1:7], spwid=[1:3], nchan=[64,64,64], start=[1,1,1], step=[1,1,1]); imgr.setimage(nx=256, ny=256, cellx='1arcsec', celly='1arcsec', stokes='I', mode='channel', nchan=30, start=47, step=4, fieldid=1, spwid=[1:3]); imgr.weight(type='natural', mosaic=T); imgr.setvp(dovp=T, usedefaultvp=T, dosquint=F); imgr.setoptions(padding=1.0, ftmachine='mosaic') imgr.setmfcontrol(scaletype="SAULT" , minpb=0.1, constpb=0.5); imgr.clean(algorithm="mfclark" , niter=300, gain=0.1, model="src.cln.model" , mask='', image="src.cln.cm" , residual="src.cln.resid", interactive=F, npercycle=30); dv.gui(); - emission in channels 5-26 - RMS in chan 4 = 39 mJy/bm chan 3 = 37 mJy/bm chan 29= 37 mJy/bm - all nice and consistent. # now clean down to 1 sigma with mask regions using mask template. imgr.clean(algorithm="mfclark" , niter=1000, gain=0.1, threshold='37mJy', model="src.cln.model2" , mask='src.cln.mask2', image="src.cln.cm2" , residual="src.cln.resid2", interactive=T, npercycle=30, masktemplate='src.cln.cm'); imgr.done(); # for channels with no emission, selected a tiny box off to the side # for channels with emission, made a different convolution region for # each channel. # Beam: 8.19528 by 5.60394 (arcsec) at pa -171.964 (deg) # Major cycles do not seem to have any improvements between them... # Maximum residuals are staying about the same. Not sure what goes. # 6 Major cycles completed. # no Glish viewer error dv.gui(); - rms and emission regions are about the same # Image just the 16apr04 data and see how the 2nd dataset affects the image: # cp -r ../16apr98/src.16apr.ms/ . imgr:=imager(filename='src.16apr.ms'); imgr.setdata(mode='channel', fieldid=[1:7], spwid=[1:3], nchan=[64,64,64], start=[1,1,1], step=[1,1,1]); imgr.setimage(nx=256, ny=256, cellx='1arcsec', celly='1arcsec', stokes='I', mode='channel', nchan=30, start=47, step=4, fieldid=1, spwid=[1:3]); imgr.weight(type='natural', mosaic=T); imgr.setvp(dovp=T, usedefaultvp=T, dosquint=F); imgr.setoptions(padding=1.0, ftmachine='mosaic') imgr.setmfcontrol(scaletype="SAULT" , minpb=0.1, constpb=0.5); imgr.clean(algorithm="mfclark" , niter=300, gain=0.1, model="src.cln.model3" , mask='', image="src.cln.cm3" , residual="src.cln.resid3", interactive=F, npercycle=30); imgr.done(); dv.gui(); - rms in chan 1 = 39 mJy/bm - rms in chan 4 = 42 mJy/bm - rms in chan 29= 39 mJy/bm RMS slightly worse but not significantly so - peak = 1.913 Jy in channel 25 - about the same as previous image - in fact, since the calibration is so poor because of the weather and the gcal is resolved, the 2nd dataset doesn't improve the image... Leave it out? - give testers the option. ----------------------------------------- # changed syntax... dach: im:=image(infile='src.cln.cm'); im2:=im.moments(outfile='src.mom0.all', moments=0, axis=4, mask='indexin(4,[5:26])', includepix=[0.070,1000.0]); RMS ~ 7 Jy/bm-km/s psfile = 0.mom0.ps - this is a rough image since I only cleaned average channels. - looks about like the MIRIAD image. im3:=im.moments(outfile='src.mom1.all', moments=1, axis=4, mask='indexin(4,[5:26])', includepix=[0.7,1000.0]); im4:=im.moments(outfile='src.mom1.all2', moments=1, axis=4, mask='indexin(4,[5:26])', includepix=[0.14,1000.0]); im5:=im.moments(outfile='src.mom1.all3', moments=1, axis=4, mask='indexin(4,[5:26])', includepix=[0.21,1000.0]); im6:=im.moments(outfile='src.mom1.all4', moments=1, axis=4, mask='indexin(4,[5:26])', includepix=[0.24,1000.0]); im7:=im.moments(outfile='src.mom1.all5', moments=1, axis=4, mask='indexin(4,[5:26])', includepix=[0.28,1000.0]); OK, got most of the background emission down to zero. Use this one. dach, have to close out all image tools generated... im.done(); im2.done(); im3.done(); im4.done(); im5.done(); im6.done(); im7.done(); Set min = 235, max = 570 km/s rainbow 1 output in 0.mom1.ps Looks reasonable. # bug(): color bar doesn't print to ps file... appears to be related to another bug when I couldn't get the dec label to print when the colorbar was present. Increasing margins to maximum on canvas manager doesn't help. Added text to current bug report requesting fix... imgr.mem(algorithm='mfentropy', niter=10, sigma='0.1Jy', constrainflux=F, mask='', model='src.mem.model', residual='src.mem.resid', image='src.mem.cm'); imgr.done(); dv.gui(); # I can see the source but, again, the mosaic doesn't appear to be # as good as the single field. What goes? -------------------------------------------------------------------- -------------------------------------------------------------------- -- testing in daily: -- can I now specify spwid[1:4] and not have imgr crash? -- imgr:=imager(filename='src.all.ms'); imgr.setdata(mode='channel', fieldid=[1:7], spwid=[1:4], nchan=[64,64,64,64], start=[1,1,1,1], step=[1,1,1,1]); imgr.setimage(nx=256, ny=256, cellx='1arcsec', celly='1arcsec', stokes='I', mode='channel', nchan=30, start=47, step=4, fieldid=1, spwid=[1:4]); imgr.weight(type='natural', mosaic=T); imgr.setvp(dovp=T, usedefaultvp=T, dosquint=F); imgr.setoptions(padding=1.0, ftmachine='mosaic') imgr.setmfcontrol(scaletype="SAULT" , minpb=0.1, constpb=0.5); imgr.clean(algorithm="mfclark" , niter=300, gain=0.1, model="src.cln.model4" , mask='', image="src.cln.cm4" , residual="src.cln.resid4", interactive=F, npercycle=30); # crash and burn - testing details in bug.summary2 imgr:=imager(filename='src.all.ms'); imgr.setdata(mode='channel', fieldid=[1:7], spwid=[1:3], nchan=[64,64,64], start=[1,1,1], step=[1,1,1]); imgr.setimage(nx=256, ny=256, cellx='1arcsec', celly='1arcsec', stokes='I', mode='channel', nchan=30, start=47, step=4, fieldid=1, spwid=[1:3]); imgr.weight(type='natural', mosaic=T); imgr.setvp(dovp=T, usedefaultvp=T, dosquint=F); imgr.setoptions(padding=1.0, ftmachine='mosaic') imgr.setmfcontrol(scaletype="SAULT" , minpb=0.1, constpb=0.5); imgr.clean(algorithm="mfclark" , niter=300, gain=0.1, model="src.cln.model6" , mask='', image="src.cln.cm6" , residual="src.cln.resid6", interactive=F, npercycle=30); dv.gui(); # crash and burn dell -r src.cln.*4 dell -r src.cln.*5 dell -r src.cln.*6 ----------------------------------------------------------------------- daily doesn't work at all - kumar fixed a bug, patched stable, now try: imgr:=imager(filename='src.all.ms'); imgr.setdata(mode='channel', fieldid=[1:7], spwid=[1:3], nchan=[64,64,64], start=[1,1,1], step=[1,1,1]); imgr.setimage(nx=256, ny=256, cellx='1arcsec', celly='1arcsec', stokes='I', mode='channel', nchan=30, start=47, step=4, fieldid=1, spwid=[1:3]); imgr.weight(type='natural', mosaic=T); imgr.setvp(dovp=T, usedefaultvp=T, dosquint=F); imgr.setoptions(padding=1.0, ftmachine='mosaic') imgr.setmfcontrol(scaletype="SAULT" , minpb=0.1, constpb=0.5); imgr.clean(algorithm="mfclark" , niter=300, gain=0.1, model="src.cln.model4" , mask='', image="src.cln.cm4" , residual="src.cln.resid4", interactive=F, npercycle=30); # OK, this is working. # try something a bit more challenging: imgr:=imager(filename='src.all.ms'); imgr.setdata(mode='channel', fieldid=[1:7], spwid=[1:4], nchan=[64,64,64,64], start=[1,1,1,1], step=[1,1,1,1]); imgr.setimage(nx=256, ny=256, cellx='1arcsec', celly='1arcsec', stokes='I', mode='channel', nchan=30, start=47, step=4, fieldid=1, spwid=[1:4]); imgr.weight(type='natural', mosaic=T); imgr.setvp(dovp=T, usedefaultvp=T, dosquint=F); imgr.setoptions(padding=1.0, ftmachine='mosaic') imgr.setmfcontrol(scaletype="SAULT" , minpb=0.1, constpb=0.5); imgr.clean(algorithm="mfclark" , niter=300, gain=0.1, model="src.cln.model4" , mask='', image="src.cln.cm4" , residual="src.cln.resid4", interactive=F, npercycle=30); # WARN: Visibility channels in spw 4 is not being used -- GOOD # humm there are 13 or 14 of these messages - we only need one. dv.gui(); --------------------------------------- Test calibrator async parm (George doesn't know what it should do and says it hasn't been tested at all). Before I put it in the cookbook - I want to know what goes with it. cal:=calibrater(filename='all.ms'); cal.reset(); cal.setdata(msselect='FIELD_ID IN [1,2] && SPECTRAL_WINDOW_ID IN [1]', async=T); cal.setsolve(type='G', t=0, refant=4, table='tstt.gcal'); cal.solve(); # WHOA: named arg "async" does not match any formal in call to function (val mode = none, val nchan = 1, val start = Looks like it doesn't work at all... #################################################################### #################################################################### #################################################################### CLEAN UP: dell -r 1230+123.ms dell -r 3c273.* dell -r ngc4826.* dell -r gcal.* dell -r *.resid* dell -r *.model* dell -r FTScreen dell -r Screen dell -r *flags* dell -r src.cln.cm2 dell -r src.cln.cm3 dell -r dell -r *~