# From - Tue May 14 11:15:29 2002 # # Hi gang, # # Appended below is the spectral line script we discussed last week. We've # sorted out the (unrelated) problems in weekly (the new one is in pretty # good shape w.r.t. this script). In Socorro, stable has been patched with # the B calibration fixes. # # Keep in mind that this is intended just as a guide as you get the hang of # it. We'll pass along some updates as defects are fixed, etc., over the # next few weeks. Stay tuned for info on how to use autoflag and how to do # the continuum subtraction. In the tutorials, you can use the final script # we converge upon as you see fit. Steve suggested a website with links # into the relevant parts of the docs; I don't think this is necessary, but # we can certainly provide a hand-out with some guidance on where to look in # the docs for basic info. # # Cheers, # George # # # 2002 Summer School Spectral Line Data Reduction Tutorial # (gmoellen 02May08) # # # Observation: VLA D-configuration, 1413 MHz, 1 IF, 63 channels, RR/LL # Fields: Target: NGC5921 (HI) (FieldId=3) # Calibrator: 1445+099 (FieldId=2) # Flux Density Calibrator: 1331+305 (FieldId1) # Standard Calibration/Imaging Sequence: # ------------------------------------- # 1. Import Data from UVFITS # 2. Summary/Examination # 3. (Editing) # 4. Set Flux Density model for 1331+305 (3C286) # 5. Solve for G (1331+305 & 1445+099) (with channel selection) # 6. Solve for B (1331+305) (applying G) # 7. Scale 1445+099 G amps according to 1331+305 G amps # 8. Form corrected data for imaging, applying B and G # 9. Form channel image (no explicit restore needed!) # 10. (Continuum subtraction) # 11. (Image analysis: profiles, fits, moments) # # Possible variations: # ------------------- # - Use B solution to get better G solutions (better F.D. scaling) # - Get B from 1445+099 data instead of 1331+305 data # - Vary channel selection for G solutions # - Vary channel selection/combination for imaging # - Vary imaging weighting options (natural, robust) # - Interactive imaging # - Self-cal on continuum source (~84mJy) in NGC5921 field (non-line chans) # - Image after each calibration step to demonstrate effect include 'imager.g'; include 'calibrater.g'; include 'image.g'; include 'ms.g'; include 'flagger.g'; include 'sysinfo.g' include 'logger.g' # # VLA H_I line data reducion: NGC 5921 # # Load the data # aipsroot:=sysinfo().root(); # Assemble name of UVFITS file fitsfilename:=spaste(aipsroot, '/data/nrao/VLA/ngc5921.fits'); msetS := fitstoms(msfile='ngc5921.ms', # Line dataset fitsfile=fitsfilename); msetS.summary(verbose=T); # Obtain summary msetS.done(); # Finish ms tool # # Basic editing (more to come, including autoflag) # fgS:=flagger('ngc5921.ms'); # Make flagger tool fgS.setflagmode('flag'); # Set mode to flag data: fgS.auto(); # Flag ACs. fgS.done(); # Finish flagger tool: # # Set source model # include 'imager.g' imgrS:=imager('ngc5921.ms'); # Start imager tool imgrS.setjy(fieldid=1, fluxdensity=[14.8009,0,0,0]); # Set 1331+305 (3C286) model # # Solve for G Jones # calS := calibrater('ngc5921.ms'); # Create calibrater tool calS.setdata(msselect='FIELD_ID <= 2', # Select data for calibrators mode='channel', # Use only ~center 40 chans start=11, nchan=40); calS.setsolve(type='G', # Arrange to solve for G on t=300.0, # 300-second timescale refant=15, table='ngc5921.gcal'); calS.state(); # Review setsolve settings calS.solve(); # Solve calS.plotcal(tablename='ngc5921.gcal'); # Inspect solutions # # Solve for B Jones (bandpass) # calS.setdata(msselect='FIELD_ID==1'); # Select BP calibrator (1331+305) calS.reset(); # Reset apply/solve state calS.setapply(type='G', # Arrange to apply G solutions t=0.0, table='ngc5921.gcal'); calS.setsolve(type='B', # Arrange to solve for bandpass t=86400.0, # (one solution for whole) refant=15, table='ngc5921.bcal'); calS.state(); # Review setapply/setsolve settings calS.solve(); # Solve calS.plotcal(tablename='ngc5921.bcal'); # Inspect # JEH: NEXT TWO LINES HAVE BEEN ADDED BY ME: # BASENAME:='ngc5921'; # include 'specline_plotbp.g'; # # Establish the flux density scale # calS.fluxscale(tablein='ngc5921.gcal', # Transfer flux density scale tableout='ngc5921.fluxcal', # from 3c286 to others reference='1331+30500002', transfer=['1445+09900002']); calS.plotcal(tablename='ngc5921.fluxcal'); # Inspect # # Correct the line data for G and B Jones # # Select fields to which calibration # will be applied calS.setdata(msselect='FIELD_ID IN [1:3]'); calS.reset(); # Reset setapply/setsolve calS.setapply(type='G', # Arrange to apply G solutions t=0.0, table='ngc5921.fluxcal', select='FIELD_ID==2'); # only solns from 1445+099 calS.setapply(type='B', # Arrange to apply B solutions t=0.0, table='ngc5921.bcal'); calS.state(); # Review setapply settings calS.correct(); # Apply solutions and write # CORRECTED_DATA column calS.done(); # Finish calibrater tool # JEH: NEXT TWO LINES HAVE BEEN ADDED BY ME: include 'msplot.g'; msplot('ngc5921.ms',edit=T); #!!! THIS BOMBS!!! Bugged, but no response. # Make a channel map of NGC 5921 # # imgrS.setdata(fieldid=3, # Select data for field 3, spectral # spwid=1, # spectral window 1, # mode='channel', # nchan=56, # avoid edge channels # start=6, # step=1); # imgrS.setimage(nx=256, # Imaging parameters # ny=256, # cellx='10arcsec', # celly='10arcsec', # stokes='I', # mode='channel', # nchan=27, # average ch in pairs # start=6, # step=2, # fieldid=3); # imgrS.weight(type='uniform'); # Uniform weighting # JEH: THOSE ARE NOT THE PARAMETERS I WOULD USE. # JEH: I WOULD USE THE FOLLOWING: imgrS.setdata(fieldid=3, # Select data for field 3, spectral spwid=1, # spectral window 1, mode='channel', nchan=63, # avoid edge channels start=1, step=1); imgrS.setimage(nx=256, # Imaging parameters ny=256, cellx='10arcsec', celly='10arcsec', stokes='I', mode='channel', nchan=63, start=1, step=1, fieldid=3); imagr.weight(type='briggs', rmode='norm', robust=+0); imgrS.clean(algorithm='clark', # Image and deconvolve niter=3000, # with Clark CLEAN threshold='0.0Jy', model='ngc5921.model', image='ngc5921.image', residual='ngc5921.residual'); imgrS.done(); # Finish imager tool # # View image in viewer: # myim:=image('ngc5921.image') myim.view() ######################################## # New stuff added by JEH: ######################################## # SUBTRACTING CONTINUUM: # ######################################## # Store the image coordinate system so that regions # can be made in world coordinates cs:=myim.coordsys(); # Define region of image which is continuum in ICHANSEL. # Note that a union of regions is a compound region that # has to be made from two *world* regions, not pixel regions, # so we have to construct it from regions made with the # wrange function from the default regionmanager (drm). # See the reference manual on regionmanager for more on this. ICHANSEL := [2,6,50,59]; ## line-free channels ## start1, stop1, start2, stop2 # Make first region cont1 for lower continuum channels range1:=dq.quantity([ICHANSEL[1],ICHANSEL[2]],'pix'); cont1:=drm.wrange(pixelaxis=4,range=range1,csys=cs); ## NOTE!! Velocity is along the 4th axis! # Make second region cont2 for upper continuum channels range2:=dq.quantity([ICHANSEL[3],ICHANSEL[4]],'pix'); cont2:=drm.wrange(pixelaxis=4,range=range2,csys=cs); ## NOTE!! Velocity is along the 4th axis! # Make region cont3 that specifies the union of the upper and # lower continuum channels cont3:=drm.union(cont1,cont2); # Form integrated continuum as moment=-1 of the line cube # only using the continuum channels myim.moments(outfile='ngc5921.cont', moments=-1, axis=4, region=cont3); # Subtract continuum from line cube, writing continuum subtracted # cube to disk: nocont:=imagecalc(outfile='ngc5921.nocont', pixels='$myim - ngc5921.cont'); # View continuum subtracted datacube nocont.view(); #################################################### # Moment analysis: ## The following replicates behaviour of MOMNT in AIPS ## ## NOTE! there are a number of very flexible options ## in the image.moment tool, and the user would do well ## to read up on it rather than treating it like the ## black-box that MOMNT in AIPS is!! ## ## (You could also try using the image.momentsgui tool ## instead - I haven't tested this). ## We will use the ICHANSEL parameter used above to ## derive our channels of line emission. These are ## used to define a region to be used in the moment analysis line_range:=dq.quantity([ICHANSEL[2]+1,ICHANSEL[3]-1],'pix'); line_region:=drm.wrange(pixelaxis=4,range=line_range,csys=cs); myim.moments(moments=[0, 1, 2] , smoothaxes=[1,2,4] , # Pixels to be included in # moment analysis will be # decided by examing the # flux in a smoothed version # of the data. # We will smooth in ra, dec # and velocity (axis=1,2,4) region = line_region, # Only use channels with # line emission mask='$myim>0', # Only include pixels with # positive flux in original map excludepix=[0.003], # Only sum flux from pixels # in smoothed map with flux # greater than abs(0.001 Jy) smoothtypes="gauss gauss hann" ,# Gaussian smoothing in ra # and dec; Hanning smoothing # in velocity smoothwidths=[5,5,3] , # Width of smoothing kernals doppler="OPTICAL" , # optical velocity definition outfile='ngc5921.nocont'); # Take a spectra through to check. # There are two interesting ways to pick positions for taking # spectra. The first is to view the continuum map and to take # spectra at locations of your continuum sources to make sure # that continuum subtraction has gone ok. The second is to # look at your moment0 map and use this to look at regions with # line emission. The latter is the default behavior, however the # entire cube has crap at the beginig and end channels, so we # explicitly use the moment0 map we just calculated. # Click on the "bullseye" icon with a mouse button to use it # for selecting the region for taking a spectra; keep it pressed to # continuously update spectra. See reference manual for more. include 'imageprofilefitter.g'; myprof1 := imageprofilefitter( infile='ngc5921.nocont', infile2='ngc5921.cont'); myprof2:=imageprofilefitter( infile='ngc5921.nocont', infile2='ngc5921.nocont.integrated'); myim.done();