Testing - NGC 7538 regression script DSS - 23feb07 ----------------------------------------------------------------------- Summary of issues: Most issues associated with the NGC 7538 regression script itself are in the notes of 22feb07. This testing deals with mosaic imaging on the 7538 calibrated data produced yesterday. Unfortunately, I was not able to get to the GG Tau dataset. - The fact that the inputs are global makes it very difficult to run scripts they way I do (e.g., cut an paste from a text file). I constantly have to use the restore() function or go into the command line to see what the inputs are. Global parameters are, in my opinion, only useful when working from the command line when you want to minimize typing. When working from a script, they are a hindrance rather than a help. Indeed, I typically loose quite a bit of time (10-20%) because I often forget to run the restore() command before each task run and then I get crazy results that just don't make sense. It can take me a while to figure out what happened. I suppose that after I get burned on this often enough, I'll learn to use the restore function before every single event in a script. However, given that this is such a problem for me, perhaps it will be for others too. Is it possible to turn off the global parameters? So everything returns to the default after each task run? Say, run autorestore('on') once to turn off the stickiness of the parameters, then run autorestore('off') to allow global parameters to remain set after a task. - In viewer, I set the scale for pixel coordinates but the axes still show RA and Dec. They should show pixel labels and coordinates. Cursor position is given in pixels however so something is working, just not the display. - During mosaicing, the logger says: NORMAL imager::weight(): Briggs weighting: sidelobes will be suppressed over full image : using 400 pixels in the uv plane This looks like the fields are not being weighted separately as they should (either that or the logger info is incorrect). If the logger message is what I think it is, CASA appears to weight all data from a single phase center (taken from the first input field). This is not the correct way to set the image weights in mosaic data - each field should have its own weighting function (which is what aips++ does when you weight the data with mosaic=T) The other disconcerting message is: NORMAL SkyEquation::fixImageScale: Using No image plane weighting Is my briggs weighting really not being applied or is this refering to something else. I think the logger message needs to be corrected if there is, in fact, image plane weighting. - I first created a mosaic image in aips++ that I was semi-reasonable Even with different deconvolution regions on each plane, the amount of cleaning per channel needed to be different, so I got a semi-good image but not something I would be willing to publish. For publication quality imaging, I'll have to process each image plane separately. But at least I got something out that was about 90% what I could reasonably expect. Then I tried to reproduce the mosaic in CASA. In the end, I could produce a mosaiced image of the N7538 dataset but it was not appropriate for this data. I needed a wider field for each image (odd that CASA cuts the mosaic beam response off more than AIPS++ and I can't control this in the mosaic task). I tried to use the tool kit to shift the phase center and see if I could manage to bring the continuum image into the mosaic a bit more but I could never get the tool kit to work. Hopefully this is user error because the documentation is not well defined for the tool kit. In any event, I have to stop testing so I don't have time to proceed further. I did find out that the clean boxes cannot be defined for each channel separately so I'll have to clean each image plane separately and define boxes for each one to get out even the first image. This will take much more work than aips++. However, I was not able to get to this point. ----------------------------------------------------------------------- ----------------------------------------------------------------------- ----------------------------------------------------------------------- Actual detailed test notes: alias casadaily 'source /home/ballista/casa/daily/casainit.csh; casapy' casadaily Fields: 6 ID Name Right Ascension Declination Epoch 0 1328+307 13:31:08.29 +30.30.33.04 J2000 1 2229+695 22:30:36.48 +69.46.28.00 J2000 2 NGC7538C 23:14:02.48 +61.27.14.86 J2000 3 NGC7538D 23:13:43.82 +61.27.00.18 J2000 4 NGC7538E 23:13:34.64 +61.27.26.44 J2000 5 NGC7538F 23:13:35.76 +61.28.33.66 J2000 Spectral Windows:(2 unique spectral windows and 2 unique polarization setups) SpwID #Chans Frame Ch1(MHz) Resoln(kHz) TotBW(kHz) Ref(MHz) Corrs 0 63 TOPO 23691.4682 118.164062 6152.34375 23694.4955 RR 1 63 TOPO 23719.6063 118.164062 6152.34375 23722.6336 LL Useful commands to work around problems: ------------------------------------------ tb.clearlocks() mp.clearplot restore() --------------- From yesterday: restore() split('ngc7538.ms','split.n7538.ms',fieldid=[2,3,4,5],spwid=[0,1], nchan=62,start=0,step=1,datacolumn='CORRECTED_DATA') # spwid = 0 should be the NH3 line at 23.6914682 GHz # primary beam at K band is 2' = 120" # synthesized beam is about 2.8" Mosaic the four fields, first figure out how to do this in app and see if I can reproduce the results in CASA: stable AIPS++: imgr:=imager(filename='split.n7538.ms'); imgr.setdata(mode='none', fieldid=[1:4],spwid=1); imgr.setimage(nx=400, ny=400, cellx='1arcsec', celly='1arcsec', stokes='I', mode='channel', nchan=55, start=3, step=1, spwid=[1],fieldid=2); imgr.weight(type='briggs',rmode='norm',robust=0.5, mosaic=T); imgr.setvp(dovp=T, usedefaultvp=T, dosquint=F); imgr.setoptions(padding=1.0, ftmachine='mosaic') imgr.clean(algorithm='hogbom', niter=1000, gain=0.1, threshold='', interactive=T, npercycle=100, model='n7538-app1.mod', residual='n7538-app1.resid', image='n7538-app1.cm', masktemplate=''); # created large masks for each channel as appropriate. # beam: 2.86171 by 2.74022 (arcsec) at pa 70.4054 (deg) dv.gui() # the emission is more extensive than I thought, I'll have to make larger mask regions. imgr.clean(algorithm='hogbom', niter=1000, gain=0.1, threshold='', interactive=T, npercycle=100, model='n7538-app2.mod', residual='n7538-app2.resid', image='n7538-app2.cm', masktemplate='n7538-app1.cm'); # rms in a plane with no line emission (just the continuum source): about 4.3 mJy/beam # I'm not cleaning deep enough. I still need to do some work to optimize the mask. The structure is so extended, multi-resolution clean would be better suited. The fields are not, by any stretch of the imagination, Nyquist sampled - they don't even have regular spacing. I could make the mosaic better centered on the fields by choosing the phase center and decreasing ny, increase nx slightly. Also, the continuum source is way off the side of one of the fields, Note, setoptions padding=1.0 = default (so I don't need to set it explicitely). imgr.done(); After quite a bit of testing, making mistakes, having crashes, summary: #imgr.setmfcontrol(scaletype="SAULT", minpb=0.05, constpb=0.4, # fluxscale='n7538-app3.fluxscale'); # SAULT weighting not good for fields that don't overlap well. # Use NONE. # there seems to be a bug in the gridding when nx does not equal ny. The image doesn't turn out the same - I have some source emission that is too close to what looks like the edge of the mosaic when nx not= ny but it is ok when nx=ny. In both cases, there are many pixels outside of the effective mosaic beam pattern so the source structure is no where close to the edge of the gridded image. imcent := dm.direction('J2000', '23h13m47.8', '61d27m40'); imgr:=imager(filename='split.n7538.ms'); imgr.setdata(mode='none', fieldid=[1:4],spwid=1); imgr.setimage(nx=400, ny=400, cellx='1arcsec', celly='1arcsec', stokes='I', mode='channel', nchan=55, start=3, step=1, spwid=[1],fieldid=2, doshift=T, phasecenter=imcent); imgr.weight(type='briggs',rmode='norm',robust=0.5, mosaic=T); imgr.setvp(dovp=T, usedefaultvp=T, dosquint=F); imgr.setoptions(ftmachine='mosaic') imgr.clean(algorithm='mfclark', niter=2000, gain=0.1, threshold='', interactive=T, npercycle=100, model='n7538-app3.mod', residual='n7538-app3.resid', mask='n7538-app3.mask', image='n7538-app3.cm', masktemplate=''); # When cleaning, after each 100 iterations, it looks like the dirty image is being displayed on the interactive clean GUI. So I can't tell how deep to clean - or maybe so little flux is being cleaned out that it looks still like the dirty image. After 200 iterations, cut my losses and let it clean the remaining 1800 iterations to see what I've got. Cleaning is not near deep enough for the continuum and some of the channels still, mask looks about right. One problem is that only 45 iterations are being done on channels with only continuum emission while many more are done on the line emission channels. Classic problem that each channel needs to have different clean parameters # clean deeper: imgr.clean(algorithm='mfclark', niter=5000, gain=0.1, threshold='', interactive=F, npercycle=100, model='n7538-app4.mod', residual='n7538-app4.resid', mask='n7538-app3.mask', image='n7538-app4.cm', masktemplate=''); imgr.done(); # SEVERE: its.tabs[ddname] is not an agent, imager crashed after the first major cycle... Try again, maybe it didn't like the fact that the mask was being displayed on the viewer while it was being used? Worked the 2nd time so maybe that was it. # so, one channel reached the max of 5000 iterations, channels with continuum only still only got to a level of 45 iterations so not enough flux was cleaned out of them. This bit about processing all channels in the mosaic with the same parameters is really foolish for this kind of image. To do better I would have to process each channel separately - I don't have time for this. # cut my losses again, stop working with AIPS++ after 5 hours and see if I can reproduce anything close to this in CASA. ------------------------------------------------------------ CASA: # Try the mosaic task - there is no documentation in the cookbook on # the mosaic task (only tool-based help). It appears that: # I can't shift the phase center... # I see no way of specifying which field should be the phase center if # I don't actually shift the phase center to a specific value. # I can only use clean boxes, in-line help says: cleanbox -- List of [blc-x,blc-y,trc-x,trc-y] values default: []; example: cleanbox=[110,110,150,145] Note: This can also be a filename with clean values: fieldindex blc-x blc-y trc-x trc-y # I have no idea how to set up clean boxes that are different for each # channel. I talked to Mark, this is not possible in AIPS and # probably not for CASA either. Try to get pixel locations for a cleanbox region. In viewer, I set the scale for pixel coordinates but the axes still show RA and Dec. The cursor position is shown in pixels so this will do for now. I have to get a first image before I can figure out cleanbox inputs - clean without boxes: restore() mosaic(vis='split.n7538.ms',imagename='n7538-casa1', mode='channel',mfalg='mfclark', niter=1000,gain=0.1,threshold=0.,mask='', nchan=55,start=3,step=1, fieldid=[0,1,2,3],spwid=[0], imsize=[400,400],cell=[1.,1.], weighting='briggs',rmode='norm',robust=0.5) # beam: 2.81312 by 2.65496 (arcsec) at pa 69.0031 (deg) # odd, this beam is slightly different than the AIPS++ beam. What goes? Are the fields being properly weighted as they should be in a mosaic or are they being weighted from a single phase center? # The logger says: Fri Feb 23 22:14:57 2007 NORMAL imager::weight(): Briggs weighting: sidelobes will be suppressed over full image : using 400 pixels in the uv plane So it looks like the fields are not being weighted separately as they should (either that or the logger info is incorrect). This won't work - apparently the first field is considered the phase center. See if I can re-order the field inputs to modify the field center to be closer to where I want it. restore() mosaic(vis='split.n7538.ms',imagename='n7538-casa2', mode='channel',mfalg='mfclark', niter=1000,gain=0.1,threshold=0.,mask='', nchan=55,start=3,step=1, fieldid=[1,2,3,0],spwid=[0], imsize=[400,400],cell=[1.,1.], weighting='briggs',rmode='norm',robust=0.5) OK, this worked, but I'll need a larger imsize, the first mosaic field is off the map. Also, the continuum source is not being cleaned - it is near the edge of the mosaic and apparently this mosaic algorithm cut off the mosaic right in the middle of the continuum source. Also, the field is being seriously overcleaned in the inner part of the mosaic - each channel should not be cleaned in the same manner. This is odd because I think I'm setting everything exactly as I did in aips++ when I did not do a phase shift. Apparently, CASA cuts off the edge of the mosaic before aips++ does. try using minpb < 0.1 which is the default. # from aips++ documentation: minpb: This is to defined up to what level the voltage pattern is going to applied when using setvp. The default is 0.1 of the primary beam or the voltage pattern defined for the antenna. So this should make the actual mosaic-ed field larger although the image itself will be the same (because there are blank pixels around the mosaic). restore() mosaic(vis='split.n7538.ms',imagename='n7538-casa3', mode='channel',mfalg='mfclark', niter=1000,gain=0.1,threshold=0.,mask='', nchan=55,start=3,step=1, fieldid=[1,2,3,0],spwid=[0], imsize=[400,400],cell=[1.,1.], weighting='briggs',rmode='norm',robust=0.5, minpb=0.02) This did not help, I can't get access to the continuum source so I won't be able to clean it. OK, it is clear that I can't get a mosaic image that works for this data. Try going to aips++ in python: # I can't figure out how to use the direction measures to define the phase center of a mosaic. There is very little documentation that I can find. The help on im says there is a phase center parameter but there is no help on it. The cookbook says it can be a direction string or field identifier. No help on how to create a "direction string." restore() imcent=dm.direction('J2000', '23h13m47.8', '61d27m40') # this won't work, dm is not defined. import dm # no module named dm... But I found a measures module, me, that has a direction function. imcent=me.direction('J2000', '23h13m47.8', '61d27m40') typing imcent gives: {'m0': {'unit': 'rad', 'value': -0.20160007301577768}, 'm1': {'unit': 'rad', 'value': 1.0726987508229593}, 'refer': 'J2000', 'type': 'direction'} # odd that 23h converts to a negative value radian. Something is wrong. And it doesn't work anyway. # Kumar said to try: imcent=me.direction('J2000, 23h13m47.8 61d27m40') typing imcent gives: {'m0': {'unit': 'rad', 'value': 0.0}, 'm1': {'unit': 'rad', 'value': 1.5707963267948966}, 'refer': 'J2000', 'type': 'direction'} - now only one of the values is given, m1 and I still get a Quantity error. Kumar said I can't use setdata and setimage, I have to use selectvis and defineimage which have different input formats and parameters... imcent=me.direction( 'J2000, 23h13m47.8, 61d27m40') im.open('split.n7538.ms'); im.selectvis(field=[0,1,2,3],spwid=0); im.defineimage(nx=400, ny=400, cellx='1arcsec', celly='1arcsec', stokes='I', mode='channel', nchan=55, start=3, step=1, spwid=0, phasecenter=imcent) im.weight(type='briggs',rmode='norm',robust=0.5, mosaic=true); im.setvp(dovp=T, usedefaultvp=T, dosquint=F); im.setoptions(ftmachine='mosaic') im.clean(algorithm='mfclark', niter=2000, gain=0.1, threshold='', model='n7538-casa4.mod', residual='n7538-casa4.resid', image='n7538-casa4.cm') : Fri Feb 23 23:51:57 2007 SEVERE MosaicFT::findConvFunction: Convolution function is misbehaved - support seems to be zero # I'm not asking for a mask, it should just clean the inner quarter of the mosaic. # maybe I have a user error in the inputs but help im.clean indicates that my inputs and their format are most likely correct (but then again, the in-line help also says I can do this with interactive clean using a masktemplate. I suspect that the in-line help is just a copy of aips++ here and it is not all implemented). # try using another input to threshold (no in-line help on this) im.clean(algorithm='mfclark', niter=2000, gain=0.1, threshold=0., model='n7538-casa4.mod', residual='n7538-casa4.resid', image='n7538-casa4.cm') Fri Feb 23 23:55:59 2007 SEVERE imager::removeTable() (file /home/ballista/casa/daily/code/synthesis/implement/MeasurementEquations/Imager.cc, line 6688): Table n7538-casa4.resid is already open in the process. It needs to be closed first Fri Feb 23 23:55:59 2007 SEVERE imager::clone() (file /home/ballista/casa/daily/code/synthesis/implement/MeasurementEquations/Imager.cc, line 5415): Exception: Invalid Table operation: SetupNewTable /home/sola/dss/casa.testing/casa/regression/ngc7538/separate.processing/n7538-casa4.resid is already opened (is in the table cache) Out[10]: False Apparently the first error put imager in a bad state so the second run of clean crashed big time. Exit and restart casapy daily: ----------------------------------------------------------- imcent=me.direction( 'J2000, 23h13m47.8, 61d27m40') im.open('split.n7538.ms'); im.selectvis(field=[0,1,2,3],spwid=0); im.defineimage(nx=400, ny=400, cellx='1arcsec', celly='1arcsec', stokes='I', mode='channel', nchan=55, start=3, step=1, spwid=0, phasecenter=imcent) im.weight(type='briggs',rmode='norm',robust=0.5, mosaic=true); im.setvp(dovp=T, usedefaultvp=T, dosquint=F); im.setoptions(ftmachine='mosaic') im.clean(algorithm='mfclark', niter=2000, gain=0.1, threshold=0., model='n7538-casa4.mod', residual='n7538-casa4.resid', image='n7538-casa4.cm') : Fri Feb 23 23:57:59 2007 SEVERE MosaicFT::findConvFunction: Convolution function is misbehaved - support seems to be zero ### I quit, I managed to get all the inputs accepted for each of the mosaic setup parameters but that is all. I can't seem to deconvolve this mosaic image. Maybe mosaic deconvolution will work for another mosaic image that is more standard. Could not produce a mosaic image with CASA that is appropriate for this source. If I get this to work, these are the box regions I'll need to input to CASA - however, I'm sure that they will make a poor image because I need vastly different deconvolution regions for each channel. The only way I might be able to do this is clean each channel with different clean box regions. But this will give me a start to begin the processing. This is from the aips++ generated image that is 400 x 400 pixels, no phase shift, alas. For continuum emission: blc = 173 263 (x,y) trc = 201 293 There is line emission in channel 27 over most of the mosaic, it shifts around extensively from channel to channel. Try: blc 1 = 142 142 2 = 40 175 3 = 227 271 trc 1 = 249 273 2 = 105 241 3 = 292 345 Most channels have no emission in these regions but I think most of the line emission is in these 3 boxes. ----------------------- ---------------------------------------------- ---------------------------------------------------------------------