casa: regressions/ngc4826-mms_tutorial_regression.py Source File

Go to the documentation of this file.
00001 ##########################################################################
00002 #                                                                        #
00003 # Demo Script for NGC 4826 (BIMA line data)                              #
00004 #                                                                        #
00005 # (with additions for Multi-MS testing, DP 2012)                         # 
00006 #                                                                        #
00007 # Converted by  STM 2008-05-27 (Beta Patch 2.0) new tasking/clean/cal    #
00008 # Updated by     CB 2008-05-30                  start from raw data      #
00009 # Updated by    STM 2008-06-01                  scriptmode, plotting     #
00010 # Updated by CB,STM 2008-06-02                  improvements             #
00011 # Updated by CB,STM 2008-06-03                  bigger cube, pbcor       #
00012 # Updated by CB,STM 2008-06-04                  pbcor stuff              #
00013 # Updated by CB,STM 2008-06-04                  tutorial script          #
00014 # Updated by CB     2008-06-05                  small tweaks             #
00015 # Updated by CB,STM 2008-06-12                  final revisions (DS)     #
00016 # Updated by STM    2008-06-14                  post-school fix          #
00017 # Updated by STM    2008-06-17                  scriptmode for regress   #
00018 # Updated by STM    2008-06-19                  regression dictionaries  #
00019 # Updated by STM    2008-06-30                  add channel check        #
00020 # Updated by STM    2008-07-08                  bug fixes                #
00021 # Updated by STM    2008-10-06                  Patch 3 regression       #
00022 # Updated by STM    2008-10-16                  Patch 3 immoments        #
00023 # Updated by STM    2008-10-22                  testing version          #
00024 # Updated by STM    2008-11-03                  minor updates            #
00025 # Updated by STM    2008-11-19                  Patch 3 400x400 image    #
00026 # Updated by STM    2008-12-01                  Patch 3 release          #
00027 # Updated by STM    2008-06-03                  Patch 4 vals, fix error  #
00028 # Updated by STM    2009-06-19                  Patch 4 release          #
00029 # Updated by STM    2009-12-07                  Release 3.0 final        #
00030 #                                                                        #
00031 # N4826 - BIMA SONG Data                                                 #
00032 #                                                                        #
00033 # This data is from the BIMA Survey of Nearby Galaxies (BIMA SONG)       #
00034 # Helfer, Thornley, Regan, et al., 2003, ApJS, 145, 259                  #
00035 # Many thanks to Michele Thornley for providing the data and description #
00036 #                                                                        #
00037 # First day of observations only                                         #
00038 # Trial for summer school                                                #
00039 # NOTE: REGRESSION VERSION FOR CASA TESTING                              #
00040 #                                                                        #
00041 # Script Notes:                                                          #
00042 #    o This script has some interactive commands, such as with plotxy    #
00043 #      and the viewer.  If scriptmode=True, then this script will stop   #
00044 #      and require a carriage-return to continue at these points.        #
00045 #    o Sometimes cut-and-paste of a series of lines from this script     #
00046 #      into the casapy terminal will get garbled (usually a single       #
00047 #      dropped character). In this case, try fewer lines, like groups    #
00048 #      of 4-6.                                                           #
00049 #    o The results are written out as a dictionary in a pickle file      #
00050 #         out.ngc4826.tutorial.regression.<datestring>.pickle            #
00051 #      as well as in a text file                                         #
00052 #         out.ngc4826.tutorial.<datestring>.log                          #
00053 #      (these are not auto-deleted at start of script)                   #
00054 #    o This script keeps internal regression values, but you can provide #
00055 #      a file ngc4826_tutorial_regression.pickle from a previous run     #
00056 #                                                                        #
00057 ##########################################################################
00058 import time
00059 import os
00060 import pickle
00061 
00062 # Crate multi-MS at which step?
00063 usemms_step    = -1 # not at all
00064 #usemms_step   = 0 # at concat
00065 #usemms_step   = 1 # before flagging
00066 #usemms_step   = 2 # before setjy
00067 #usemms_step   = 3 # before gaincal
00068 #usemms_step   = 4 # before fluxscale
00069 #usemms_step   = 5 # before applycal
00070 #usemms_step   = 6 # before split
00071 #usemms_step   = 7 # before clean
00072 
00073 
00074 # Use field-wise or spw-wise splitting in multi-MS?
00075 usefieldwise = True
00076 
00077 
00078 def makefieldwisemms(msname):    
00079     casalog.post('******* Making field-wise MMS. *******', "INFO", 'makefieldwisemms')
00080     print '******* Making field-wise MMS. *******'
00081     os.system('rm -rf temp.ms; mv '+msname+' temp.ms')
00082     members = []
00083     tb.open('temp.ms/FIELD')
00084     nfield = tb.nrows()
00085     tb.close()
00086     for i in range(0,nfield):
00087         member = 'ngc4826.tutorial-field'+str(i)+'.ms'
00088         os.system('rm -rf '+member)
00089         print member
00090         split(vis='temp.ms', outputvis=member, field=i, datacolumn='all')
00091         members.append(member)
00092     concat(vis=members, concatvis=msname, createmms=True)
00093 
00094 def makespwwisemms(msname):
00095     casalog.post('******* Making spw-wise MMS. *******', "INFO", 'makespwwisemms')
00096     print '******* Making spw-wise MMS. *******'
00097     os.system('rm -rf temp.ms; mv '+msname+' temp.ms')
00098     members = []
00099     tb.open('temp.ms/SPECTRAL_WINDOW')
00100     nspw = tb.nrows()
00101     tb.close()
00102     for i in range(0,nspw):
00103         member = 'ngc4826.tutorial-spw'+str(i)+'.ms'
00104         os.system('rm -rf '+member)
00105         print member
00106         split(vis='temp.ms', outputvis=member, spw=str(i), datacolumn='all')
00107         members.append(member)
00108     concat(vis=members, concatvis=msname, createmms=True)
00109 
00110 
00111 #scriptmode = True
00112 
00113 # If you want to run like a regression, including making PNG for plots,
00114 # then set to False
00115 scriptmode = False
00116 
00117 # Enable benchmarking?
00118 benchmarking = True
00119 
00120 # Sets a shorthand for fixed input script/regression files
00121 scriptprefix='ngc4826_tutorial_regression'
00122 
00123 #
00124 ##########################################################################
00125 #                                                                        
00126 # Clear out previous run results
00127 os.system('rm -rf ngc4826.tutorial.*')
00128 
00129 # Sets a shorthand for the ms, not necessary
00130 prefix='ngc4826.tutorial'
00131 msfile = prefix + '.16apr98.ms'
00132 
00133 print 'Tutorial Regression Script for BIMASONG NGC4826 Mosaic'
00134 print 'Version for Release 0 (3.0.0) 7-Dec-2009'
00135 print 'Will do: import, flagging, calibration, imaging'
00136 print ''
00137 #
00138 ##########################################################################
00139 #
00140 # 
00141 ##########################################################################
00142 #
00143 # N4826 - BIMA SONG Data
00144 # 16apr98
00145 #       source=ngc4826
00146 #       phasecal=1310+323
00147 #       fluxcal=3c273, Flux = 23 Jy on 16apr98
00148 #       passcal= none - data were observed with online bandpass correction.
00149 #
00150 # NOTE: This data has been filled into MIRIAD, line-length correction 
00151 #       done, and then exported as separate files for each source.
00152 #       3c273 was not line length corrected since it was observed
00153 #       for such a short amount of time that it did not need it.  
00154 #
00155 # From miriad: source Vlsr = 408; delta V is 20 km/s 
00156 #
00157 #
00158 ##########################################################################
00159 #
00160 if benchmarking:
00161     startTime=time.time()
00162     startProc=time.clock()
00163 
00164 _mydatapath = os.environ.get('CASAPATH').split()[0]+'/data/regression/ngc4826/'
00165 
00166 ##########################################################################
00167 # Import and concatenate sources
00168 ##########################################################################
00169 #
00170 # USB spectral windows written separately by miriad for 16apr98
00171 # Assumes these are in sub-directory called "fitsfiles" of working directory
00172 print '--Importuvfits (16apr98)--'
00173 default('importuvfits')
00174 
00175 print "Starting from the uvfits files exported by miriad"
00176 print "The USB spectral windows were written separately by miriad for 16apr98"
00177 
00178 importuvfits(fitsfile=_mydatapath+'fitsfiles/3c273.fits5', vis='ngc4826.tutorial.3c273.5.ms')
00179 
00180 importuvfits(fitsfile=_mydatapath+'fitsfiles/3c273.fits6', vis='ngc4826.tutorial.3c273.6.ms')
00181 
00182 importuvfits(fitsfile=_mydatapath+'fitsfiles/3c273.fits7', vis='ngc4826.tutorial.3c273.7.ms')
00183 
00184 importuvfits(fitsfile=_mydatapath+'fitsfiles/3c273.fits8', vis='ngc4826.tutorial.3c273.8.ms')
00185 
00186 importuvfits(fitsfile=_mydatapath+'fitsfiles/1310+323.ll.fits9', vis='ngc4826.tutorial.1310+323.ll.9.ms')
00187 
00188 importuvfits(fitsfile=_mydatapath+'fitsfiles/1310+323.ll.fits10', vis='ngc4826.tutorial.1310+323.ll.10.ms')
00189 
00190 importuvfits(fitsfile=_mydatapath+'fitsfiles/1310+323.ll.fits11', vis='ngc4826.tutorial.1310+323.ll.11.ms')
00191 
00192 importuvfits(fitsfile=_mydatapath+'fitsfiles/1310+323.ll.fits12', vis='ngc4826.tutorial.1310+323.ll.12.ms')
00193 
00194 importuvfits(fitsfile=_mydatapath+'fitsfiles/1310+323.ll.fits13', vis='ngc4826.tutorial.1310+323.ll.13.ms')
00195 
00196 importuvfits(fitsfile=_mydatapath+'fitsfiles/1310+323.ll.fits14', vis='ngc4826.tutorial.1310+323.ll.14.ms')
00197 
00198 importuvfits(fitsfile=_mydatapath+'fitsfiles/1310+323.ll.fits15', vis='ngc4826.tutorial.1310+323.ll.15.ms')
00199 
00200 importuvfits(fitsfile=_mydatapath+'fitsfiles/1310+323.ll.fits16', vis='ngc4826.tutorial.1310+323.ll.16.ms')
00201 
00202 importuvfits(fitsfile=_mydatapath+'fitsfiles/ngc4826.ll.fits5', vis='ngc4826.tutorial.ngc4826.ll.5.ms')
00203 
00204 importuvfits(fitsfile=_mydatapath+'fitsfiles/ngc4826.ll.fits6', vis='ngc4826.tutorial.ngc4826.ll.6.ms')
00205 
00206 importuvfits(fitsfile=_mydatapath+'fitsfiles/ngc4826.ll.fits7', vis='ngc4826.tutorial.ngc4826.ll.7.ms')
00207 
00208 importuvfits(fitsfile=_mydatapath+'fitsfiles/ngc4826.ll.fits8', vis='ngc4826.tutorial.ngc4826.ll.8.ms')
00209 
00210 if benchmarking:
00211     import2time=time.time()
00212 
00213 #
00214 ##########################################################################
00215 #
00216 print '--Concat--'
00217 default('concat')
00218 
00219 concat(vis=['ngc4826.tutorial.3c273.5.ms',
00220             'ngc4826.tutorial.3c273.6.ms',
00221             'ngc4826.tutorial.3c273.7.ms',
00222             'ngc4826.tutorial.3c273.8.ms',
00223             'ngc4826.tutorial.1310+323.ll.9.ms',
00224             'ngc4826.tutorial.1310+323.ll.10.ms',
00225             'ngc4826.tutorial.1310+323.ll.11.ms',
00226             'ngc4826.tutorial.1310+323.ll.12.ms',
00227             'ngc4826.tutorial.1310+323.ll.13.ms',
00228             'ngc4826.tutorial.1310+323.ll.14.ms',
00229             'ngc4826.tutorial.1310+323.ll.15.ms',
00230             'ngc4826.tutorial.1310+323.ll.16.ms',
00231             'ngc4826.tutorial.ngc4826.ll.5.ms',
00232             'ngc4826.tutorial.ngc4826.ll.6.ms',
00233             'ngc4826.tutorial.ngc4826.ll.7.ms',
00234             'ngc4826.tutorial.ngc4826.ll.8.ms'],
00235        concatvis='ngc4826.tutorial.ms',
00236        freqtol="",dirtol="1arcsec",async=False, createmms=(usemms_step==0))
00237 
00238 if benchmarking:
00239     concat2time=time.time()
00240 
00241 #
00242 ##########################################################################
00243 #
00244 # TUTORIAL NOTES:
00245 #
00246 # You can invoke tasks in two ways:
00247 #
00248 # (1) As function calls with arguments as shown above for concat and used
00249 #     extensively in this script, e.g.
00250 #
00251 #        task( par1=val1, par2=val2, ... )
00252 #
00253 #     with parameters set as arguments in the call.  Note that in this
00254 #     case, the global parameter values are NOT used or changed, and any
00255 #     task parameters that are not specified as arguments to the call
00256 #     will be defaulted to the task-specific default values (see the
00257 #     "help task" description).
00258 #
00259 # (2) By setting the values of the global parameters and then using the
00260 #     "go" command (if taskname is set) or calling the task with no
00261 #     arguments.  For example:
00262 #
00263 #        default task
00264 #        par1 = val1
00265 #        par2 = val2
00266 #        ...
00267 #        inp
00268 #        task()
00269 #
00270 #     In this case, the "default" command sets the parmeters to their
00271 #     task defaults, and sets the "taskname" paramter to the task to be
00272 #     run.  The "inp" command displays the current values for the task
00273 #     parameters.  Then the call with no arguments runs with the globals.
00274 #
00275 #     Warning: "go" does not work inside scripts. See Cookbook.
00276 #
00277 # Using the concat call above as an example, we would do:
00278 #
00279 #default('concat')
00280 #
00281 #vis = ['ngc4826.tutorial.3c273.5.ms',
00282 #       'ngc4826.tutorial.3c273.6.ms',
00283 #       'ngc4826.tutorial.3c273.7.ms',
00284 #       'ngc4826.tutorial.3c273.8.ms',
00285 #       'ngc4826.tutorial.1310+323.ll.9.ms',
00286 #       'ngc4826.tutorial.1310+323.ll.10.ms',
00287 #       'ngc4826.tutorial.1310+323.ll.11.ms',
00288 #       'ngc4826.tutorial.1310+323.ll.12.ms',
00289 #       'ngc4826.tutorial.1310+323.ll.13.ms',
00290 #       'ngc4826.tutorial.1310+323.ll.14.ms',
00291 #       'ngc4826.tutorial.1310+323.ll.15.ms',
00292 #       'ngc4826.tutorial.1310+323.ll.16.ms',       
00293 #       'ngc4826.tutorial.ngc4826.ll.5.ms',
00294 #       'ngc4826.tutorial.ngc4826.ll.6.ms',
00295 #       'ngc4826.tutorial.ngc4826.ll.7.ms',
00296 #       'ngc4826.tutorial.ngc4826.ll.8.ms']
00297 #
00298 #concatvis='ngc4826.tutorial.ms'
00299 #freqtol = ""
00300 #dirtol = "1arcsec"
00301 #async=False
00302 #
00303 #concat()
00304 
00305 #
00306 ##########################################################################
00307 #
00308 # Fix up the MS (temporary, changes to importfits underway)
00309 #
00310 #print '--Fixing up spw rest frequencies in MS--'
00311 #vis='ngc4826.tutorial.ms'
00312 #tb.open(vis+'/SOURCE',nomodify=false)
00313 #spwid=tb.getcol('SPECTRAL_WINDOW_ID')
00314 #spwid.setfield(-1,int)
00315 # Had to do this for 64bit systems 08-Jul-2008
00316 ##spwid.setfield(-1,'int32')
00317 ##tb.putcol('SPECTRAL_WINDOW_ID',spwid)
00318 ##tb.close()
00319 
00320 # This ensures that the rest freq will be found for all spws. 
00321 
00322 #
00323 ##########################################################################
00324 # 16 APR Calibration
00325 ##########################################################################
00326 print '--Clearcal--'
00327 print 'Create scratch columns and initialize in '+'ngc4826.tutorial.ms'
00328 
00329 # Force create/initialize of scratch columns
00330 # NOTE: plotxy will not run properly without this step.
00331 #
00332 clearcal(vis='ngc4826.tutorial.ms', addmodel=False)
00333 
00334 if benchmarking:
00335     clearcal2time=time.time()
00336 
00337 #
00338 ##########################################################################
00339 #
00340 # List contents of MS
00341 #
00342 print '--Listobs--'
00343 listobs(vis='ngc4826.tutorial.ms')
00344 
00345 # Should see the listing in logger
00346 # Some parts shown below
00347 #
00348 
00349 print "There are 3 fields observed in a total of 16 spectral windows"
00350 print "   field=0    3c273    spwids 0,1,2,3               64 chans "
00351 print "   field=1    1310+323 spwids 4,5,6,7,8,9,10,11     32 chans "
00352 print "   field=2~8  NGC4826  spwids 12,13,14,15           64 chans "
00353 print ""
00354 print "See listobs summary in logger"
00355 
00356 if benchmarking:
00357     list2time=time.time()
00358 
00359 # Fields: 9
00360 # ID   Code Name          Right Ascension  Declination   Epoch   
00361 # 0         3C273         12:29:06.70      +02.03.08.60  J2000   
00362 # 1         1310+323      13:10:28.66      +32.20.43.78  J2000   
00363 # 2         NGC4826       12:56:44.24      +21.41.05.10  J2000   
00364 # 3         NGC4826       12:56:41.08      +21.41.05.10  J2000   
00365 # 4         NGC4826       12:56:42.66      +21.41.43.20  J2000   
00366 # 5         NGC4826       12:56:45.82      +21.41.43.20  J2000   
00367 # 6         NGC4826       12:56:47.39      +21.41.05.10  J2000   
00368 # 7         NGC4826       12:56:45.82      +21.40.27.00  J2000   
00369 # 8         NGC4826       12:56:42.66      +21.40.27.00  J2000   
00370 # Spectral Windows:
00371 # SpwID  #Chans Frame Ch1(MHz)    ChanWid(kHz)TotBW(kHz)  Ref(MHz)    Corrs
00372 # 0          64 LSRK  115108.478  1562.5      100000      115108.478  YY  
00373 # 1          64 LSRK  115198.615  1562.5      100000      115198.615  YY  
00374 # 2          64 LSRK  115288.478  1562.5      100000      115288.478  YY  
00375 # 3          64 LSRK  115378.615  1562.5      100000      115378.615  YY  
00376 # 4          32 LSRK  114974.256  3125        100000      114974.256  YY  
00377 # 5          32 LSRK  115074.393  3125        100000      115074.393  YY  
00378 # 6          32 LSRK  115174.256  3125        100000      115174.256  YY  
00379 # 7          32 LSRK  115274.393  3125        100000      115274.393  YY  
00380 # 8          32 LSRK  115374.256  3125        100000      115374.256  YY  
00381 # 9          32 LSRK  115474.392  3125        100000      115474.392  YY  
00382 # 10         32 LSRK  115574.255  3125        100000      115574.255  YY  
00383 # 11         32 LSRK  115674.392  3125        100000      115674.392  YY  
00384 # 12         64 LSRK  114950.191  1562.5      100000      114950.191  YY  
00385 # 13         64 LSRK  115040.205  1562.5      100000      115040.205  YY  
00386 # 14         64 LSRK  115129.946  1562.5      100000      115129.946  YY  
00387 # 15         64 LSRK  115219.96   1562.5      100000      115219.96   YY  
00388 
00389 #
00390 ##########################################################################
00391 # Plotting and Flagging
00392 ##########################################################################
00393 #
00394 # The plotxy task is the interactive x-y display and flagging GUI
00395 #
00396 #print '--Plotxy--'
00397 #default(plotxy)
00398 
00399 # First look at amplitude as a funciton of uv-distance using an
00400 # average over all times and channels for each source.
00401 #if scriptmode:
00402 #    plotxy(vis='ngc4826.tutorial.ms',xaxis='uvdist',yaxis='amp',
00403 #           field='0',spw='0~3',
00404 #           averagemode='vector',timebin='1e7',width='1000',crossscans=True,
00405 #           selectplot=True,title='Field 0 SPW 0~3')
00406 
00407 #    print "Looking at 3C273 versus uvdist with time and chan average"
00408     # Pause script if you are running in scriptmode
00409 #    user_check=raw_input('Return to continue script\n')
00410 #else:
00411 #    plotxy(vis='ngc4826.tutorial.ms',xaxis='uvdist',yaxis='amp',
00412 #           field='0',spw='0~3',
00413 #           averagemode='vector',timebin='1e7',width='1000',crossscans=True,
00414 #           selectplot=True,title='Field 0 SPW 0~3',
00415 #           interactive=False,
00416 #           figfile='ngc4826.tutorial.ms.plotxy.field0.ampuv.allavg.png')
00417 
00418 # NOTE: width here needs to be larger than combination of all channels
00419 # selected with spw and/or field. Since field and spw are unique in this
00420 # case, both don't need to be specified, however plotting is much faster
00421 # if you "help it" by selecting both.
00422 
00423 # Now lets look at the target source, the first of the NGC4826 mosaic fields
00424 # which are 2~8 in this MS.
00425 #
00426 # Since we are plotting versus velocity we can clearly see the bad edge
00427 # channels and the overlap between spw
00428 #
00429 # There is nothing terribly wrong with this data and again we will flag the
00430 # edge channels non-interactively later for consistency.
00431 #
00432 # Normally, if there were obviously bad data, you would flag it here
00433 # before calibration.  To do this, hit the Mark Region button, then draw a box around
00434 # some of the moderately high outliers, and then Flag.
00435 #
00436 # But this data is relatively clean, and flagging will not improve results.
00437 #
00438 # Interactive plotxy
00439 #if scriptmode:
00440 #    plotxy(vis='ngc4826.tutorial.ms',xaxis='velocity',yaxis='amp',
00441 #           field='2',spw='12~15',
00442 #           averagemode='vector',timebin='1e7',crossscans=True,
00443 #           selectplot=True,newplot=False,title='Field 2 SPW 12~15')
00444 
00445 #    print "You could Mark Region around outliers and Flag"
00446     # Pause script if you are running in scriptmode
00447 #    user_check=raw_input('Return to continue script\n')
00448 #else:
00449     # Output to file
00450     # Set up a Python loop to do all the N4826 fields:
00451 #    for fld in range(2,9):
00452 #        field = str(fld)
00453 #        plotxy(vis='ngc4826.tutorial.ms',xaxis='velocity',yaxis='amp',
00454 #               field=field,spw='12~15',
00455 #               averagemode='vector',timebin='1e7',crossscans=True,
00456 #               selectplot=True,newplot=False,title='Field 2 SPW 12~15',
00457 #               interactive=False,
00458 #               figfile='ngc4826.tutorial.ms.plotxy.field2.ampvel.tavg.png')
00459 
00460     # Now the 1310+323 field
00461 #    plotxy(vis='ngc4826.tutorial.ms',xaxis='velocity',yaxis='amp',
00462 #           field='1',spw='4~11',
00463 #           averagemode='vector',timebin='1e7',crossscans=True,
00464 #           selectplot=True,newplot=False,title='Field 1 SPW 4~11',
00465 #           interactive=False,
00466 #           figfile='ngc4826.tutorial.ms.plotxy.field1.ampvel.tavg.png')
00467 
00468     # Now the 3C273 field
00469     # This one should be time and channel averaged to test this
00470 #   plotxy(vis='ngc4826.tutorial.ms',xaxis='velocity',yaxis='amp',
00471 #         field='0',spw='0~3',
00472 #         averagemode='vector',timebin='1e7',crossscans=True,
00473 #           selectplot=True,newplot=False,title='Field 0 SPW 0~3',
00474 #           interactive=False,
00475 #           figfile='ngc4826.tutorial.ms.plotxy.field0.ampvel.tavg.png')
00476 
00477 # You can also have it iterate over baselines, using Next to advance
00478 # Set to NOT plot autocorrelations
00479 # Example using 3C273: (interactive only)
00480 #if scriptmode:
00481 #    plotxy(vis='ngc4826.tutorial.ms',xaxis='channel',yaxis='amp',
00482 #           field='0',spw='0~3',
00483 #           selectdata=True,antenna='*&*',
00484 #           averagemode='vector',timebin='1e7',crossscans=True,
00485 #           iteration='baseline',
00486 #           selectplot=True,newplot=False,title='Field 0 SPW 0~3')
00487         
00488     # Pause script if you are running in scriptmode
00489 #    user_check=raw_input('Return to continue script\n')
00490 
00491 #
00492 # Finally, look for bad data. Here we look at field 8 w/o averaging
00493 #if scriptmode:
00494 #    plotxy(vis='ngc4826.tutorial.ms',xaxis='time',yaxis='amp',field='8',spw='12~15',
00495 #           selectplot=True,newplot=False,title='Field 8 SPW 12~15')
00496 
00497 #    print "You can see some bad data here"
00498 #    print "Mark Region and Locate, look in logger"
00499 #    print "This is a correlator glitch in baseline 3-9 at 06:19:30"
00500 #    print "PLEASE DON\'T FLAG ANYTHING HERE. THE SCRIPT WILL DO IT!"
00501 #    print "In a normal session you could Mark Region and Flag."
00502 #    print "Here we will use flagdata instead."
00503     # Pause script if you are running in scriptmode
00504 #    user_check=raw_input('Return to continue script\n')
00505 
00506     # If you change xaxis='channel' you see its all channels
00507 #else:
00508     # Plot to file
00509 #    plotxy(vis='ngc4826.tutorial.ms',xaxis='time',yaxis='amp',
00510 #           field='8',spw='12~15',
00511 #           selectplot=True,newplot=False,title='Field 8 SPW 12~15',
00512 #           interactive=False,
00513 #           figfile='ngc4826.tutorial.ms.plotxy.field2.amptime.noavg.png')
00514     
00515 if benchmarking:
00516     plotxy2time=time.time()
00517 
00518 if (usemms_step==1):
00519     if usefieldwise:
00520         makefieldwisemms('ngc4826.tutorial.ms')
00521     else:
00522         makespwwisemms('ngc4826.tutorial.ms')
00523 
00524 
00525 #
00526 ##########################################################################
00527 #
00528 # Flag end channels
00529 #
00530 print '--Flagdata--'
00531 default('tflagdata')
00532 
00533 print ""
00534 print "Flagging edge channels in all spw"
00535 print "  0~3:0~1;62~63 , 4~11:0~1;30~31, 12~15:0~1;62~63 "
00536 print ""
00537 
00538 tflagdata(vis='ngc4826.tutorial.ms', mode='manual',
00539          spw='0~3:0;1;62;63,4~11:0;1;30;31,12~15:0;1;62;63')
00540 
00541 #
00542 # Flag correlator glitch
00543 #
00544 print ""
00545 print "Flagging bad correlator field 8 antenna 3&9 spw 15 all channels"
00546 print "  timerange 1998/04/16/06:19:00.0~1998/04/16/06:20:00.0"
00547 print ""
00548 
00549 tflagdata(vis='ngc4826.tutorial.ms', mode='manual', field='8', spw='15', antenna='3&9', 
00550          timerange='1998/04/16/06:19:00.0~1998/04/16/06:20:00.0')
00551 
00552 #
00553 # Flag non-fringing antenna 6
00554 #
00555 # NOTE: this is already flagged in the data so do nothing more here
00556 #flagdata(vis='ngc4826.tutorial.ms', mode='manualflag', antenna='6',
00557 #        timerange='1998/04/16/09:42:39.0~1998/04/16/10:24:46.0')
00558 
00559 #
00560 #
00561 ##########################################################################
00562 #
00563 # Some example clean-up editing
00564 # Slightly high almost-edge channel in field='1', spw='4' (channel 2)
00565 # can be flagged interactively with plotxy.
00566 #
00567 #plotxy(vis='ngc4826.tutorial.ms',
00568 #       xaxis='channel',yaxis='amp',field='1',spw='4',
00569 #       averagemode='vector',timebin='1e7',crossscans=True,
00570 #       selectplot=True,newplot=False,title='Field 1 SPW 4')
00571 
00572 #
00573 ##########################################################################
00574 #
00575 # Use Flagmanager to save a copy of the flags so far
00576 #
00577 print '--Flagmanager--'
00578 default('flagmanager')
00579 
00580 print "Now will use flagmanager to save a copy of the flags we just made"
00581 print "These are named myflags"
00582 
00583 flagmanager(vis='ngc4826.tutorial.ms',mode='save',versionname='myflags',
00584             comment='My flags',merge='replace')
00585 
00586 # Can also use Flagmanager to list all saved versions
00587 #
00588 flagmanager(vis='ngc4826.tutorial.ms',mode='list')
00589 
00590 if benchmarking:
00591     flag2time=time.time()
00592 
00593 print "Completed pre-calibration flagging"
00594 
00595 if (usemms_step==2):
00596     if usefieldwise:
00597         makefieldwisemms('ngc4826.tutorial.ms')
00598     else:
00599         makespwwisemms('ngc4826.tutorial.ms')
00600 
00601 
00602 #
00603 ##########################################################################
00604 #
00605 # CALIBRATION
00606 #
00607 ##########################################################################
00608 #
00609 # Bandpasses are very flat because of observing mode used (online bandpass
00610 # correction) so bandpass calibration is unnecessary for these data.
00611 #
00612 ##########################################################################
00613 #
00614 # Derive gain calibration solutions.
00615 # We will use VLA-like G (per-scan) calibration:
00616 #
00617 ##########################################################################
00618 #
00619 # Set the flux density of 3C273 to 23 Jy
00620 #
00621 print '--Setjy (3C273)--'
00622 default('setjy')
00623 
00624 setjy(vis='ngc4826.tutorial.ms',field='0',fluxdensity=[23.0,0.,0.,0.],spw='0~3', scalebychan=False, usescratch=False)
00625 #
00626 # Not really necessary to set spw but you get lots of warning messages if
00627 # you don't
00628 
00629 if benchmarking:
00630     setjy2time=time.time()
00631 
00632 if (usemms_step==3):
00633     if usefieldwise:
00634         makefieldwisemms('ngc4826.tutorial.ms')
00635     else:
00636         makespwwisemms('ngc4826.tutorial.ms')
00637 
00638 
00639 #
00640 ##########################################################################
00641 #
00642 # Gain calibration
00643 #
00644 print '--Gaincal--'
00645 default('gaincal')
00646 
00647 # This should be combining all spw for the two calibrators for single
00648 # scan-based solutions
00649 
00650 print 'Gain calibration for fields 0,1 and spw 0~11'
00651 print 'Using solint=inf combining over spw'
00652 print 'Output table ngc4826.tutorial.16apr98.gcal'
00653 
00654 gaincal(vis='ngc4826.tutorial.ms', caltable='ngc4826.tutorial.16apr98.gcal',
00655         field='0,1', spw='0~11', gaintype='G', minsnr=2.0,
00656         refant='ANT5', gaincurve=False, opacity=0.0,
00657         solint='inf', combine='spw')
00658 
00659 if benchmarking:
00660     gaincal2time=time.time()
00661 
00662 if (usemms_step==4):
00663     if usefieldwise:
00664         makefieldwisemms('ngc4826.tutorial.ms')
00665     else:
00666         makespwwisemms('ngc4826.tutorial.ms')
00667 
00668 
00669 #
00670 ##########################################################################
00671 #
00672 # Transfer the flux density scale:
00673 #
00674 print '--Fluxscale--'
00675 default('fluxscale')
00676 
00677 print ''
00678 print 'Transferring flux of 3C273 to sources: 1310+323'
00679 print 'Output table ngc4826.tutorial.16apr98.fcal'
00680 
00681 fluxscale(vis='ngc4826.tutorial.ms', caltable='ngc4826.tutorial.16apr98.gcal',
00682           fluxtable='ngc4826.tutorial.16apr98.fcal',
00683           reference='3C273', transfer=['1310+323'])
00684 
00685 # Flux density for 1310+323 is: 1.48 +/- 0.016 (SNR = 90.6, nAnt= 8)
00686 
00687 if benchmarking:
00688     fluxscale2time=time.time()
00689 
00690 if (usemms_step==5):
00691     if usefieldwise:
00692         makefieldwisemms('ngc4826.tutorial.ms')
00693     else:
00694         makespwwisemms('ngc4826.tutorial.ms')
00695 
00696 
00697 #
00698 ##########################################################################
00699 #
00700 # Plot calibration
00701 print '--Plotcal (fluxscale)--'
00702 default(plotcal)
00703 
00704 if scriptmode:
00705     # Interactive plotcal
00706     plotcal(caltable='ngc4826.tutorial.16apr98.fcal', yaxis='amp', field='')
00707     print ''
00708     print 'Plotting final scaled gain calibration table'
00709     print 'First amp vs. time for all fields '
00710         
00711     # Pause script if you are running in scriptmode
00712     user_check=raw_input('Return to continue script\n')
00713 
00714     plotcal(caltable='ngc4826.tutorial.16apr98.fcal', yaxis='phase', field='')
00715     print ''
00716     print 'and phase vs. time '
00717 
00718     # Pause script if you are running in scriptmode
00719     user_check=raw_input('Return to continue script\n')
00720 
00721     # And you can plot the SNR of the solution
00722     plotcal(caltable='ngc4826.tutorial.16apr98.fcal', yaxis='snr', field='')
00723 else:
00724     # You can also plotcal to file
00725     plotcal(caltable='ngc4826.tutorial.16apr98.fcal',yaxis='amp',field='',
00726             showgui=False,figfile='ngc4826.tutorial.16apr98.fcal.plotcal.amp.png')
00727     plotcal(caltable='ngc4826.tutorial.16apr98.fcal',yaxis='phase',field='',
00728             showgui=False,figfile='ngc4826.tutorial.16apr98.fcal.plotcal.phase.png')
00729     plotcal(caltable='ngc4826.tutorial.16apr98.fcal',yaxis='snr',field='',
00730             showgui=False,figfile='ngc4826.tutorial.16apr98.fcal.plotcal.snr.png')
00731 
00732 if benchmarking:
00733     plotcal2time=time.time()
00734 
00735 
00736 #
00737 ##########################################################################
00738 #
00739 # Correct the calibrater/target source data:
00740 # Use new parm spwmap to apply gain solutions derived from spwid1
00741 # to all other spwids... 
00742 print '--Applycal--'
00743 default('applycal')
00744 
00745 print 'Applying calibration table ngc4826.tutorial.16apr98.fcal to data'
00746 
00747 applycal(vis='ngc4826.tutorial.ms',
00748          field='', spw='',
00749          gaincurve=False, opacity=0.0, 
00750          gaintable='ngc4826.tutorial.16apr98.fcal',
00751          spwmap=[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0])
00752 
00753 if benchmarking:
00754     correct2time=time.time()
00755 
00756 if (usemms_step==6):
00757     if usefieldwise:
00758         makefieldwisemms('ngc4826.tutorial.ms')
00759     else:
00760         makespwwisemms('ngc4826.tutorial.ms')
00761 
00762 
00763 
00764 #
00765 ##########################################################################
00766 #
00767 # Check calibrated data
00768 #print '--Plotxy--'
00769 #default(plotxy)
00770 
00771 #
00772 # Here we plot the first of the NGC4826 fields unaveraged versus velocity
00773 # Notice how the spw fit together
00774 
00775 # Interactive plotxy
00776 #print "Will plot all the NGC4826 calibrated data unaveraged - this will take a while"
00777 #plotxy(vis='ngc4826.tutorial.ms',xaxis='velocity',yaxis='amp',field='2~8',spw='12~15',
00778 #       averagemode='vector',datacolumn='corrected',
00779 #       selectplot=True,newplot=False,title='Field 2~8 SPW 12~15')
00780 
00781 #print ""
00782 #print "Look for outliers, flag them if there are any bad ones"
00783 #print ""
00784         
00785 # Pause script if you are running in scriptmode
00786 #user_check=raw_input('Return to continue script\n')
00787 
00788 # You can also plot all the N4826 fields 2 through 8, for example using a loop:
00789 
00790 #for fld in range(2,9):
00791 #       field = str(fld)
00792 #       plotxy(vis='ngc4826.tutorial.ms',xaxis='velocity',yaxis='amp',
00793 #              field=field,spw='11~15',
00794 #              averagemode='vector',datacolumn='corrected',
00795 #              selectplot=True,newplot=False,title='Field '+field+' SPW 11~15')
00796 #       
00797 #       user_check=raw_input('Return to continue script\n')
00798 
00799 # Now here we time-average the data, plotting versus velocity
00800 
00801 #plotxy(vis='ngc4826.tutorial.ms',xaxis='velocity',yaxis='amp',field=field,spw=spw,
00802 #       averagemode='vector',datacolumn='corrected',
00803 #       timebin='1e7',crossscans=True,plotcolor='blue',
00804 #       selectplot=True,newplot=False,title='Field '+field+' SPW '+spw)
00805 #print ""
00806 #print 'Final Spectrum field '+field+' spw '+spw+' TimeAverage Corrected Data'
00807         
00808 # Pause script if you are running in scriptmode
00809 #user_check=raw_input('Return to continue script\n')
00810 
00811 # Here we overplot 3C273 the Time+Chan averaged calibrated and uncalibrated data
00812 
00813 #
00814 # First the corrected column in blue
00815 #field = '0'
00816 #spw = '0~3'
00817 #plotxy(vis='ngc4826.tutorial.ms',xaxis='uvdist',yaxis='amp',field=field,spw=spw,
00818 #       averagemode='vector',width='1000',datacolumn='corrected',
00819 #       timebin='1e7',crossscans=True,plotcolor='blue',
00820 #       selectplot=True,newplot=False,title='Field '+field+' SPW '+spw)
00821 #print ""
00822 #print 'Plotting field '+field+' spw '+spw+' TimeChanAverage Corrected Data in blue'
00823 #
00824 # Now the original data column in red
00825 #plotxy(vis='ngc4826.tutorial.ms',xaxis='uvdist',yaxis='amp',field=field,spw=spw,
00826 #       averagemode='vector',width='1000',datacolumn='data',
00827 #       timebin='1e7',crossscans=True,plotcolor='red',overplot=True,
00828 #       selectplot=True,newplot=False,title='Field '+field+' SPW '+spw)
00829 #
00830 #print 'OverPlotting field '+field+' spw '+spw+' TimeChanAverage Original Data in red'
00831         
00832 ## Pause script if you are running in scriptmode
00833 #user_check=raw_input('Return to continue script\n')
00834 
00835 #if benchmarking:
00836 #    plotfinal2time=time.time()
00837 
00838 print "Done calibration and plotting"
00839 #
00840 ##########################################################################
00841 #
00842 # SPLIT THE DATA INTO SINGLE-SOURCE MS
00843 # AND THEN IMAGE THE CALIBRATOR
00844 #
00845 ##########################################################################
00846 #
00847 #
00848 # Split out calibrated target source and calibrater data:
00849 #
00850 print '--Split--'
00851 default('split')
00852 
00853 print 'Splitting 3C273 data to ngc4826.tutorial.16apr98.3C273.split.ms'
00854 
00855 split(vis='ngc4826.tutorial.ms',
00856       outputvis='ngc4826.tutorial.16apr98.3C273.split.ms',
00857       field='0',spw='0~3:0~63', datacolumn='corrected')
00858 
00859 print 'Splitting 1310+323 data to ngc4826.tutorial.16apr98.1310+323.split.ms'
00860 
00861 split(vis='ngc4826.tutorial.ms',
00862       outputvis='ngc4826.tutorial.16apr98.1310+323.split.ms',
00863       field='1', spw='4~11:0~31', datacolumn='corrected')
00864 
00865 print 'Splitting NGC4826 data to ngc4826.tutorial.16apr98.src.split.ms'
00866 
00867 split(vis='ngc4826.tutorial.ms',
00868       outputvis='ngc4826.tutorial.16apr98.src.split.ms',
00869       field='2~8', spw='12~15:0~63',
00870       datacolumn='corrected')
00871 
00872 if benchmarking:
00873     split2time=time.time()
00874 
00875 #
00876 ##########################################################################
00877 #
00878 #print '--Clearcal (split data)--'
00879 
00880 # If you want to use plotxy before cleaning to look at the split ms
00881 # then you will need to force the creation of the scratch columns
00882 # (clean will also do this)
00883 #clearcal('ngc4826.tutorial.16apr98.1310+323.split.ms')
00884 
00885 #clearcal('ngc4826.tutorial.16apr98.src.split.ms')
00886 
00887 #Then you might look at the data
00888 #clearplot()
00889 #plotxy(vis='ngc4826.tutorial.16apr98.src.split.ms',xaxis='time',yaxis='amp')
00890 
00891 #if benchmarking:
00892 #    clearfinal2time=time.time()
00893 
00894 #
00895 ##########################################################################
00896 #
00897 # You might image the calibrater data:
00898 #
00899 #print '--Clean (1310+323)--'
00900 #default('clean')
00901 #
00902 #
00903 #clean(vis='ngc4826.tutorial.16apr98.1310+323.split.ms', 
00904 #      imagename='ngc4826.tutorial.16apr98.cal.clean',
00905 #      cell=[1.,1.],imsize=[256,256],
00906 #      field='0',spw='0~7',threshold=10.,
00907 #      mode='mfs',psfmode='clark',niter=100,stokes='I')
00908 
00909 # You can look at this in the viewer
00910 #viewer('ngc4826.tutorial.16apr98.cal.clean.image')
00911 
00912 #if benchmarking:
00913 #    cleancal2time=time.time()
00914 
00915 
00916 if (usemms_step==7):
00917     if usefieldwise:
00918         makefieldwisemms('ngc4826.tutorial.16apr98.3C273.split.ms')
00919         makefieldwisemms('ngc4826.tutorial.16apr98.1310+323.split.ms')
00920         makefieldwisemms('ngc4826.tutorial.16apr98.src.split.ms')
00921     else:
00922         makespwwisemms('ngc4826.tutorial.16apr98.3C273.split.ms')
00923         makespwwisemms('ngc4826.tutorial.16apr98.1310+323.split.ms')
00924         makespwwisemms('ngc4826.tutorial.16apr98.src.split.ms')
00925 
00926 
00927 #
00928 #
00929 ##########################################################################
00930 #
00931 # IMAGING OF NGC4826 MOSAIC
00932 #
00933 ##########################################################################
00934 #
00935 #          Mosaic field spacing looks like:
00936 #
00937 #          F3 (field 3)         F2 (field 2)
00938 #
00939 #   F4 (field 4)      F0 (field 0)        F1 (field 1)
00940 #
00941 #          F5 (field 5)         F6 (field 6)
00942 #
00943 # 4x64 channels = 256 channels 
00944 #
00945 # Primary Beam should be about 1.6' FWHM (7m dishes, 2.7mm wavelength)
00946 # Resolution should be about 5-8"
00947 ##########################################################################
00948 #
00949 # Image the target source mosaic:
00950 #
00951 print '--Clean (NGC4826)--'
00952 default('clean')
00953 
00954 # Make image big enough so mosaic is in inner quarter (400x400)
00955 clnsize = 400
00956 print " Creating CLEAN image of size "+str(clnsize)
00957 
00958 clean(vis='ngc4826.tutorial.16apr98.src.split.ms',
00959       imagename='ngc4826.tutorial.16apr98.src.clean',
00960       field='0~6',spw='0~3',
00961       cell=[1.,1.],imsize=[clnsize,clnsize],
00962       stokes='I',
00963       mode='channel',nchan=36,start=34,width=4,
00964       interpolation='nearest',
00965       psfmode='clark',imagermode='mosaic',ftmachine='mosaic',
00966       scaletype='SAULT',
00967 ### As we moved to clean by default in flat sigma rather than
00968 ### flat snr it converges less well
00969  ###     cyclefactor=1.5,niter=10000,threshold='45mJy',
00970       cyclefactor=4,niter=10000,threshold='45mJy',
00971       minpb=0.3,pbcor=False, usescratch=False)
00972 
00973 ### NOTE: mosaic data ...Sault weighting implies a noise unform image
00974 
00975 ### Using ftmachine='mosaic', can also use ftmachine='ft' for more
00976 ### traditional image plane mosaicing
00977 
00978 ### NOTE: that niter is set to large number so that stopping point is
00979 ### controlled by threshold.
00980 
00981 ### NOTE: with pbcor=False, the final image is not "flux correct",
00982 ### instead the image has constant noise despite roll off in power as
00983 ### you move out from the phase center(s). Though this format makes it
00984 ### "look nicest", for all flux density measurements, and to get an
00985 ### accurate integrated intensity image, one needs to divide the
00986 ### srcimage.image/srcimage.flux in order to correct for the mosaic
00987 ### response pattern. One could also achieve this by setting pbcor=True
00988 ### in clean.
00989 
00990 # Try running clean adding the parameter interactive=True.
00991 # This parameter will periodically bring up the viewer to allow
00992 # interactive clean boxing. For poor uv-coverage, deep negative bowls
00993 # from missing short spacings, this can be very important to get correct
00994 # integrated flux densities.
00995 
00996 if benchmarking:
00997     clean2time=time.time()
00998 
00999 #
01000 ##########################################################################
01001 #
01002 # Do interactive viewing of clean image
01003 print '--Viewer--'
01004 if scriptmode:
01005     viewer('ngc4826.tutorial.16apr98.src.clean.image')
01006 
01007     print ""
01008     print "This is the non-pbcorrected cube of NGC4826"
01009     print "Use tape deck to move through channels"
01010     print "Close the viewer when done"
01011     print ""
01012 
01013     # Pause script if you are running in scriptmode
01014     user_check=raw_input('Return to continue script\n')
01015 
01016 #
01017 ##########################################################################
01018 #
01019 # Statistics on clean image cube
01020 #
01021 print '--ImStat (Clean cube)--'
01022 
01023 srcstat = imstat('ngc4826.tutorial.16apr98.src.clean.image')
01024 
01025 print "Found image max = "+str(srcstat['max'][0])
01026 
01027 # 256x256: refpix = '128,128,128,128'
01028 # 400x400: refpix = '200,200,200,200'
01029 # 512x512: refpix = '256,256,256,256'
01030 # 800x800: refpix = '400,400,400,400'
01031 refpix = int(clnsize/2)
01032 refbox = str(refpix)+','+str(refpix)+','+str(refpix)+','+str(refpix)
01033 print "  Using Reference Pixel "+refbox
01034 
01035 # 256x256: offbox = '106,161,153,200'
01036 # 400x400: offbox = '178,233,225,272'
01037 # 512x512: offbox = '234,289,281,328'
01038 # 800x800: offbox = '378,433,425,472'
01039 blcx = refpix - 22
01040 blcy = refpix + 33
01041 trcx = refpix + 25
01042 trcy = refpix + 72
01043 offbox = str(blcx)+','+str(blcy)+','+str(trcx)+','+str(trcy)
01044 print "  Using Off-Source Box "+offbox
01045 
01046 offstat = imstat('ngc4826.tutorial.16apr98.src.clean.image',
01047                  box=offbox)
01048 
01049 print "Found off-source image rms = "+str(offstat['sigma'][0])
01050 
01051 # 256x256: cenbox = '108,108,148,148'
01052 # 400x400: cenbox = '180,180,320,320'
01053 # 512x512: cenbox = '236,236,276,276'
01054 # 800x800: cenbox = '380,380,420,420'
01055 blcx = refpix - 20
01056 blcy = refpix - 20
01057 trcx = refpix + 20
01058 trcy = refpix + 20
01059 cenbox = str(blcx)+','+str(blcy)+','+str(trcx)+','+str(trcy)
01060 print "  Using On-Source Box "+cenbox
01061 
01062 # offlinechan = '0,1,2,3,4,5,30,31,32,33,34,35'
01063 
01064 offlinestat = imstat('ngc4826.tutorial.16apr98.src.clean.image',
01065                      box=cenbox,
01066                      chans='0,1,2,3,4,5,30,31,32,33,34,35')
01067 
01068 print "Found off-line image rms = "+str(offlinestat['sigma'][0])
01069 
01070 #
01071 ##########################################################################
01072 #
01073 # Statistics on clean model
01074 #
01075 print '--ImStat (Clean model)--'
01076 
01077 modstat = imstat('ngc4826.tutorial.16apr98.src.clean.model')
01078 
01079 print "Found total model flux = "+str(modstat['sum'][0])
01080 
01081 #
01082 ##########################################################################
01083 #
01084 # Manually correct for mosaic response pattern using .image/.flux images
01085 #
01086 print '--ImMath (PBcor)--'
01087 
01088 immath(outfile='ngc4826.tutorial.16apr98.src.clean.pbcor',
01089        mode='evalexpr',
01090        expr="'ngc4826.tutorial.16apr98.src.clean.image'/'ngc4826.tutorial.16apr98.src.clean.flux'")
01091 
01092 # now pbcor the model, be careful to mask zeros
01093 
01094 immath(outfile='ngc4826.tutorial.16apr98.src.clean.pbcormod',
01095        mode='evalexpr',
01096        expr="'ngc4826.tutorial.16apr98.src.clean.model'['ngc4826.tutorial.16apr98.src.clean.model'!=0.0]/'ngc4826.tutorial.16apr98.src.clean.flux'")
01097 
01098 
01099 #
01100 ##########################################################################
01101 #
01102 # Statistics on PBcor image cube
01103 #
01104 print '--ImStat (PBcor cube)--'
01105 
01106 pbcorstat = imstat('ngc4826.tutorial.16apr98.src.clean.pbcor')
01107 
01108 print "Found image max = "+str(pbcorstat['max'][0])
01109 
01110 pbcoroffstat = imstat('ngc4826.tutorial.16apr98.src.clean.pbcor',
01111                       box=offbox)
01112 
01113 print "Found off-source image rms = "+str(pbcoroffstat['sigma'][0])
01114 
01115 pbcorofflinestat = imstat('ngc4826.tutorial.16apr98.src.clean.pbcor',
01116                           box=cenbox,
01117                           chans='0,1,2,3,4,5,30,31,32,33,34,35')
01118 
01119 print "Found off-line image rms = "+str(pbcorofflinestat['sigma'][0])
01120 
01121 #
01122 # Statistics on PBcor image cube
01123 #
01124 print '--ImStat (PSF)--'
01125 
01126 psfstat = imstat('ngc4826.tutorial.16apr98.src.clean.psf',
01127                  box=refbox,chans='27')
01128 
01129 print "Found PSF value at refpix = "+str(psfstat['mean'][0])+" (should be 1.0)"
01130 
01131 if benchmarking:
01132     math2time=time.time()
01133 
01134 #
01135 # Statistics on PBcor model cube
01136 #
01137 print '--ImStat (PBcor model)--'
01138 
01139 pbcormodstat = imstat('ngc4826.tutorial.16apr98.src.clean.pbcormod')
01140 
01141 print "Found total model flux = "+str(pbcormodstat['sum'][0])
01142 
01143 #
01144 ##########################################################################
01145 #
01146 # Do zeroth and first moments
01147 #
01148 # NGC4826 LSR velocity is 408 km/s; delta is 20 km/s
01149 
01150 # NOTE: before 02-Jul-2008 (5631) the planes were 1-based, are now 0-based
01151 # was planes 7~28, now 6~27
01152 
01153 print '--ImMoments--'
01154 default('immoments')
01155 
01156 momfile = 'ngc4826.tutorial.16apr98.moments'
01157 momzeroimage = 'ngc4826.tutorial.16apr98.moments.integrated'
01158 momoneimage = 'ngc4826.tutorial.16apr98.moments.mom1'
01159 
01160 print "Calculating Moments 0,1 for PBcor image"
01161 
01162 # In the following we will let immoments figure out which axis
01163 # to collapse along, the spectral axis=3
01164 
01165 immoments(imagename='ngc4826.tutorial.16apr98.src.clean.pbcor',
01166           moments=[0],
01167           chans='6~27',
01168           outfile='ngc4826.tutorial.16apr98.moments.integrated') 
01169 
01170 # TUTORIAL NOTES: For moment 0 we use the image corrected for the
01171 # mosaic response to get correct integrated flux densities. However,
01172 # in *real signal* regions, the value of moment 1 is not dependent on
01173 # the flux being correct so the non-pb corrected SAULT image can be
01174 # used, this avoids having lots of junk show up at the edges of your
01175 # moment 1 image due to the primary beam correction. Try it both ways
01176 # and see for yourself.
01177 
01178 # TUTORIAL NOTES:
01179 #
01180 # Moments greater than zero need to have a conservative lower
01181 # flux cutoff to produce sensible results.
01182 
01183 immoments(imagename='ngc4826.tutorial.16apr98.src.clean.image',
01184           moments=[1],includepix=[0.2,1000.0],
01185           chans='6~27',
01186           outfile='ngc4826.tutorial.16apr98.moments.mom1') 
01187 
01188 # Now view the resulting images
01189 if scriptmode:
01190     viewer('ngc4826.tutorial.16apr98.moments.integrated')
01191     #
01192     print "Now viewing Moment-0 ngc4826.tutorial.16apr98.moments.integrated"
01193     print "Note PBCOR effects at field edge"
01194     print "Change the colorscale to get better image"
01195     print "You can also Open and overlay Contours of Moment-1 ngc4826.tutorial.16apr98.moments.mom1"
01196     print "Close the viewer when done"
01197     
01198     # Pause script if you are running in scriptmode
01199     user_check=raw_input('Return to continue script\n')
01200 
01201 # Do a moment one on channel 0 to check that the indexing is right
01202 # NOTE: THIS STILL CRASHES
01203 try:
01204     immoments(imagename='ngc4826.tutorial.16apr98.src.clean.image',
01205               moments=[1],includepix=[],
01206               chans='0',
01207               outfile='ngc4826.tutorial.16apr98.moments.plane0.mom1') 
01208 except:
01209     pass
01210 
01211 # Do a moment one on channel 35 to check that the indexing is right
01212 try:
01213     immoments(imagename='ngc4826.tutorial.16apr98.src.clean.image',
01214           moments=[1],includepix=[],
01215           chans='35',
01216           outfile='ngc4826.tutorial.16apr98.moments.plane35.mom1')
01217 except:
01218     pass
01219 
01220 if benchmarking:
01221     moments2time=time.time()
01222 
01223 #
01224 ##########################################################################
01225 #
01226 # Statistics on moment images
01227 #
01228 print '--ImStat (Moment images)--'
01229 
01230 momzerostat=imstat('ngc4826.tutorial.16apr98.moments.integrated')
01231 
01232 try:
01233     print "Found moment 0 max = "+str(momzerostat['max'][0])
01234     print "Found moment 0 rms = "+str(momzerostat['rms'][0])
01235 except:
01236     pass
01237 
01238 momonestat = imstat('ngc4826.tutorial.16apr98.moments.mom1')
01239 
01240 try:
01241     print "Found moment 1 median = "+str(momonestat['median'][0])
01242 except:
01243     pass
01244 
01245 
01246 
01247 ia.open('ngc4826.tutorial.16apr98.src.clean.image')
01248 csys=ia.coordsys()
01249 vel0=0.0
01250 vel35=0.0
01251 
01252 try:
01253     momoneplane0=imstat('ngc4826.tutorial.16apr98.moments.plane0.mom1')
01254     print "Found plane 0 moment 1 value = "+str(momoneplane0['median'][0])
01255 except:
01256     pass
01257 
01258 
01259 try:
01260     momoneplane35=imstat('ngc4826.tutorial.16apr98.moments.plane35.mom1')
01261     print "Found plane 35 moment 1 value = "+str(momoneplane35['median'][0])
01262 except:
01263     pass
01264 
01265 if(type(momoneplane0)==bool):
01266     vel0=csys.frequencytovelocity(ia.toworld([0,0,0,0])['numeric'][3])
01267 if(type(momoneplane35)==bool):
01268     vel35=csys.frequencytovelocity(ia.toworld([0,0,0,35])['numeric'][3])
01269 
01270 
01271 #
01272 ##########################################################################
01273 #
01274 # Get MS stats
01275 #
01276 ms.open('ngc4826.tutorial.16apr98.1310+323.split.ms')
01277 vismean_cal=pl.mean(ms.range(["amplitude"]).get("amplitude"))
01278 ms.close()
01279 ms.open('ngc4826.tutorial.16apr98.src.split.ms')
01280 vismean_src=pl.mean(ms.range(["amplitude"]).get("amplitude"))
01281 ms.close()
01282 
01283 #
01284 ##########################################################################
01285 #
01286 # An alternative is to mask the pbcor image before calculating
01287 # moments.  The following block shows how to do this.
01288 
01289 # To do this, open the clean pbcor file in the viewer and use the
01290 # Region Manager to create a region file
01291 
01292 #print '--Viewer--'
01293 #viewer(ngc4826.tutorial.16apr98.src.clean.pbcor)
01294 
01295 # TUTORIAL NOTES: After loading ngc4826.tutorial.16apr98.src.clean.pbcor
01296 # in the viewer as a raster, click on the file icon in top left
01297 # corner, and select the momzero image but open as a Contour map
01298 # instead of Raster. Then decrease "contour scale factor" in "Data
01299 # Display Options" gui to something like 10. select "region manager
01300 # tool" from "tool" drop down menu In region manager tool select "all
01301 # axes". Then assign the sqiggly Polygon button to a mouse button by
01302 # clicking on it with a mouse button. Then draw a polygon region
01303 # around galaxy emission, avoiding edge regions. Then in "region
01304 # manager tool" select "save last region".
01305 
01306 # Pause script if you are running in scriptmode
01307 #user_check=raw_input('Return to continue script\n')
01308 
01309 # You should have created region file ngc4826.tutorial.16apr98.src.clean.pbcor.rgn
01310 
01311 # Now use immath to use the region file to mask the cube
01312 #
01313 #print '--ImMath (masking)--'
01314 #immath(outfile='ngc4826.tutorial.16apr98.src.clean.pbcor.masked',
01315 #       mode='evalexpr',
01316 #       expr="'ngc4826.tutorial.16apr98.src.clean.pbcor'",
01317 #       region='ngc4826.tutorial.16apr98.src.clean.pbcor.rgn')
01318 
01319 #
01320 # And then make masked moment images
01321 
01322 #print '--ImMoments (masked)--'
01323 #print 'Creating masked moment 0 image ngc4826.tutorial.16apr98.moments.integratedmasked'
01324 #       
01325 #immoments(imagename='ngc4826.tutorial.16apr98.src.clean.pbcor.masked',
01326 #          moments=0,axis=3,
01327 #          planes='6~27',
01328 #          outfile='ngc4826.tutorial.16apr98.moments.integratedmasked') 
01329 #
01330 #print 'Creating masked moment 1 image ngc4826.tutorial.16apr98.moments.mom1masked'
01331 #
01332 #immoments(imagename='ngc4826.tutorial.16apr98.src.clean.pbcor.masked',
01333 #          moments=1,axis=3,
01334 #          includepix=[0.2,1000.0],
01335 #          planes='6~27',
01336 #          outfile='ngc4826.tutorial.16apr98.moments.mom1masked') 
01337 
01338 # Now view the resulting images
01339 #viewer('ngc4826.tutorial.16apr98.moments.integratedmasked')
01340 #
01341 #print "Now viewing masked Moment-0 ngc4826.tutorial.16apr98.moments.integratedmasked"
01342 #print "You can Open and overlay Contours of Moment-1 ngc4826.tutorial.16apr98.moments.mom1masked"
01343 #
01344 # Pause script if you are running in scriptmode
01345 #user_check=raw_input('Return to continue script\n')
01346 
01347 # Finally, can compute and print statistics
01348 #print '--ImStat (masked moments)--'
01349 #
01350 #maskedmomzerostat = imstat('ngc4826.tutorial.16apr98.moments.integratedmasked')
01351 #print "Found masked moment 0 max = "+str(maskedmomzerostat['max'][0])
01352 #print "Found masked moment 0 rms = "+str(maskedmomzerostat['rms'][0])
01353 #
01354 #maskedmomonestat=imstat('ngc4826.tutorial.16apr98.moments.mom1masked')
01355 #print "Found masked moment 1 median = "+str(maskedmomonestat['median'][0])
01356 
01357 if benchmarking:
01358     endProc=time.clock()
01359     endTime=time.time()
01360 
01361 #
01362 ##########################################################################
01363 # Previous results to be used for regression
01364 
01365 # Canonical regression values (using this script STM 2008-06-04) were:
01366 #vis_mean_cal=4.3269
01367 #vis_mean_src=156.992
01368 #clean_image_max=1.45868253708
01369 #clean_offsrc_rms=0.0438643493782
01370 #clean_offline_rms=0.0544108718199
01371 #clean_momentzero_max=169.420959473
01372 #clean_momentzero_rms=14.3375244141
01373 # And from STM 2008-06-18
01374 #clean_momentone_median=428.43
01375 # And from STM 2008-06-30
01376 #clean_momentone_planezero = 688.575012
01377 #clean_momentone_planelast = 119.659264
01378 # MS mean STM 2008-07-02
01379 #vis_mean_cal = 194.613642
01380 #vis_mean_src = 54.583024
01381 
01382 #New values STM 2008-10-06 Patch3 (gaincal fix)
01383 #clean_image_max = 1.466937
01384 #clean_offsrc_rms = 0.043549
01385 #clean_offline_rms = 0.054727
01386 #clean_momentzero_max = 174.315887
01387 #clean_momentzero_rms = 14.601299
01388 #clean_momentone_median = 427.726257
01389 #clean_momentone_planezero = 688.575012
01390 #clean_momentone_planelast = 119.659264
01391 #vis_mean_cal = 194.914215
01392 #vis_mean_src = 54.626986
01393 
01394 #New model fluxes STM 2008-10-23 Patch3
01395 #  from run w/o clean mosaic pb change using 256x256 image
01396 #model_sum = 73.183142
01397 #model_pbcor_sum = 76.960971
01398 
01399 #New values STM 2008-12-01 Patch3.0 (released version)
01400 #for 400x400 clean 
01401 #testdate = '2008-12-01 (STM)'
01402 #testvers = 'CASA Version 2.3.0 Rev 6654'
01403 #clean_image_max = 1.481322
01404 #clean_offsrc_rms = 0.043665
01405 #clean_offline_rms = 0.055379
01406 #clean_momentzero_max = 174.731827
01407 #clean_momentzero_rms = 14.858011
01408 #clean_momentone_median = 428.499237
01409 #clean_momentone_planezero = 688.575012
01410 #clean_momentone_planelast = 119.659264
01411 #vis_mean_cal = 194.915497
01412 #vis_mean_src = 54.627127
01413 #model_sum = 72.437971
01414 #model_pbcor_sum = 70.417830
01415 
01416 #New values STM 2009-02-25 Patch3.1
01417 #for 400x400 clean 
01418 #testdate = '2009-02-25 (STM)'
01419 #testvers = 'CASA Version 2.3.1 Rev 6826'
01420 #clean_image_max = 1.481322
01421 #clean_offsrc_rms = 0.043665
01422 #clean_offline_rms = 0.055379
01423 #clean_momentzero_max = 174.731842
01424 #clean_momentzero_rms = 14.858011
01425 #clean_momentone_median = 428.499237
01426 #clean_momentone_planezero = 688.575012
01427 #clean_momentone_planelast = 119.659264
01428 #vis_mean_cal = 194.915497
01429 #vis_mean_src = 54.627127
01430 #model_sum = 72.437973
01431 #model_pbcor_sum = 70.417831
01432 
01433 #New values STM 2009-06-19 Patch4 (released version)
01434 #for 400x400 clean
01435 #testdate = '2009-06-19 (STM)'
01436 #testvers = 'CASA Version 2.4.0 Rev 8115'
01437 #clean_image_max = 1.418478
01438 #clean_offsrc_rms = 0.043584
01439 #clean_offline_rms = 0.056824
01440 #clean_momentzero_max = 171.601685
01441 #clean_momentzero_rms = 15.532441
01442 #clean_momentone_median = 428.499237
01443 #clean_momentone_planezero = 688.575012
01444 #clean_momentone_planelast = 119.659264
01445 #vis_mean_cal = 194.915497
01446 #vis_mean_src = 54.627127
01447 #model_sum = 70.707405
01448 #model_pbcor_sum = 63.006854
01449 
01450 #New values KG 2009-11-01 Release 0 (prerelease version)
01451 #for 400x400 clean
01452 #new values for flat noise clean
01453 #testdate = '2009-12-02 (KG)'
01454 #testvers = 'CASA Version 3.0.0 Rev 9692'
01455 #clean_image_max = 1.46
01456 #clean_offsrc_rms = 0.0573
01457 #clean_offline_rms = 0.05429
01458 #clean_momentzero_max = 165.7
01459 #clean_momentzero_rms = 15.1
01460 #clean_momentone_median = 422.84
01461 #clean_momentone_planezero = 688.575012
01462 #clean_momentone_planelast = 119.659264
01463 #vis_mean_cal = 194.915497
01464 #vis_mean_src = 54.627127
01465 #model_sum = 74.374
01466 #model_pbcor_sum = 65.319
01467 
01468 ## #New values STM 2009-12-02 Release 0 (prerelease version)
01469 ## #for 400x400 clean
01470 ## #new values for flat noise clean
01471 ## testdate = '2009-12-02 (STM)'
01472 ## testvers = 'CASA Version 3.0.1 Rev 10130'
01473 ## clean_image_max = 1.465047 
01474 ## clean_offsrc_rms = 0.058497
01475 ## clean_offline_rms = 0.055416
01476 ## clean_momentzero_max = 163.726852
01477 ## clean_momentzero_rms = 15.206372
01478 ## clean_momentone_median = 428.326385
01479 ## clean_momentone_planezero = 696.702393
01480 ## clean_momentone_planelast = 127.786629
01481 ## vis_mean_cal = 194.915085
01482 ## vis_mean_src = 54.627020
01483 ## model_sum = 71.171693
01484 ## model_pbcor_sum = 61.853749
01485 
01486 ## #New values RR 2010-03-30 3.0.1 prerelease
01487 ## #for 400x400 clean
01488 ## #new values after start channel change.  I did not update passing values.
01489 ## testdate = '2010-03-30 (RR)'
01490 ## testvers = 'CASA Version 3.0.1 (build #10841)'
01491 ## clean_image_max = 1.615747 # was 1.465047 Peak hits a channel better?
01492 ## clean_offsrc_rms = 0.058497
01493 ## clean_offline_rms = 0.055416
01494 ## clean_momentzero_max = 163.726852
01495 ## clean_momentzero_rms = 15.206372
01496 ## clean_momentone_median = 429.658844 # was 428.326385; change << 1 chanwidth.
01497 ## clean_momentone_planezero = 696.702393
01498 ## clean_momentone_planelast = 127.786629
01499 ## vis_mean_cal = 194.915085
01500 ## vis_mean_src = 54.627020
01501 ## model_sum = 71.171693
01502 ## model_pbcor_sum = 66.882499 # was 61.853749 Peak hits a channel better?
01503 
01504 ## # slight change in regression values - reason not known yet.
01505 ## # Remy's changes to cleanhelper, or something from Sanjay or Kumar?
01506 ## testdate = '2010-04-24 (RR)'
01507 ## testvers = 'CASA Version 3.0.2 (build #11181)'
01508 ## clean_image_max = 1.615747
01509 ## clean_offsrc_rms = 0.058497
01510 ## clean_offline_rms = 0.055416
01511 ## clean_momentzero_max = 163.726852
01512 ## clean_momentzero_rms = 15.206372
01513 ## clean_momentone_median = 423.6954 # was 429.6588; change << 1 chanwidth.
01514 ## clean_momentone_planezero = 696.702393
01515 ## clean_momentone_planelast = 127.786629
01516 ## vis_mean_cal = 194.915085
01517 ## vis_mean_src = 54.627020
01518 ## model_sum = 71.171693
01519 ## model_pbcor_sum = 75.92 # was 66.88 Peak hits a channel better?
01520 
01521 # slight change in 1 regression value - Kumar fixed a bug in setjy with
01522 # multiple spws, which could affect this script.
01523 testdate = '2010-04-29 (RR)'
01524 testvers = 'CASA Version 3.0.2 (build #11306)'
01525 clean_image_max = 1.615747
01526 clean_offsrc_rms = 0.058497
01527 clean_offline_rms = 0.055416
01528 clean_momentzero_max = 163.726852
01529 clean_momentzero_rms = 15.206372
01530 #
01531 #  32 bits gets 423.6954 and 64 bits gets 422.142792
01532 #  diff << 1 chanwidth.
01533 clean_momentone_median = 422.92
01534 clean_momentone_planezero = 696.702393
01535 clean_momentone_planelast = 127.786629
01536 vis_mean_cal = 195.0509
01537 vis_mean_src = 54.665
01538 model_sum = 75.7516
01539 model_pbcor_sum = 75.92 # was 66.88 Peak hits a channel better?
01540 
01541 # RR, 1/19/2011 - The rmses went down, just like in orionmos4sim.  This is
01542 # good, so I won't complain too loudly.  The moment 1 median and pbcor_sum have
01543 # jiggled around a fair bit.  The median is _not_ affected by the two spurious
01544 # blobs at 501.64 km/s, though.  (Verified by doing imstat with a tight polygon
01545 # region.)
01546 clean_offsrc_rms = 0.04736
01547 clean_offline_rms = 0.049904
01548 clean_momentzero_rms = 13.737679
01549 # The chanwidth is ~16 km/s.
01550 clean_momentone_median = 435.368103
01551 
01552 model_pbcor_sum = 77.7232
01553 
01554 ## # RR, 3/11/2011 - The rmses went up, but not to their historical maxima.  The
01555 ## # model_pbcor_sum went down, but not to its historical minimum.  Nobody seems
01556 ## # to know why.
01557 ## clean_offsrc_rms = 0.0535
01558 ## clean_offline_rms = 0.0563
01559 ## model_pbcor_sum = 62.5022
01560 ## # RR, 3/12/2011 - And now, mysteriously, they're back to their 1/19/2011 values.
01561 
01562 ## RR, 3/25 - 4/3/2011, after clean was changed to used the center of output
01563 ## channel frequencies, instead of center of the first input channel in each
01564 ## output channel.
01565 clean_image_max = 1.4637
01566 clean_momentone_median = 424.3254
01567 clean_momentone_planezero = 690.6068
01568 clean_momentone_planelast = 121.6911
01569 
01570 canonical = {}
01571 canonical['exist'] = True
01572 
01573 canonical['date'] = testdate
01574 canonical['version'] = testvers
01575 canonical['user'] = 'smyers'
01576 canonical['host'] = 'sandrock'
01577 canonical['cwd'] = '/home/sandrock/smyers/Testing/Patch4/N4826'
01578 print "Using internal regression from "+canonical['version']+" on "+canonical['date']
01579 
01580 canonical_results = {}
01581 canonical_results['clean_image_max'] = {}
01582 canonical_results['clean_image_max']['value'] = clean_image_max
01583 canonical_results['clean_image_offsrc_max'] = {}
01584 canonical_results['clean_image_offsrc_max']['value'] = clean_offsrc_rms
01585 canonical_results['clean_image_offline_max'] = {}
01586 canonical_results['clean_image_offline_max']['value'] = clean_offline_rms
01587 canonical_results['clean_momentzero_max'] = {}
01588 canonical_results['clean_momentzero_max']['value'] = clean_momentzero_max
01589 canonical_results['clean_momentzero_rms'] = {}
01590 canonical_results['clean_momentzero_rms']['value'] = clean_momentzero_rms
01591 # And from STM 2008-06-18
01592 canonical_results['clean_momentone_median'] = {}
01593 canonical_results['clean_momentone_median']['value'] = clean_momentone_median
01594 # And from STM 2008-06-30
01595 canonical_results['clean_momentone_planezero'] = {}
01596 canonical_results['clean_momentone_planezero']['value'] = clean_momentone_planezero
01597 canonical_results['clean_momentone_planelast'] = {}
01598 canonical_results['clean_momentone_planelast']['value'] = clean_momentone_planelast
01599 canonical_results['clean_psfcenter'] = {}
01600 canonical_results['clean_psfcenter']['value'] = 1.0
01601 
01602 # MS mean STM 2008-07-02
01603 canonical_results['vis_mean_cal'] = {}
01604 canonical_results['vis_mean_cal']['value'] = vis_mean_cal
01605 canonical_results['vis_mean_src'] = {}
01606 canonical_results['vis_mean_src']['value'] = vis_mean_src
01607 
01608 # Model fluxes STM 2008-10-22
01609 canonical_results['model_sum'] = {}
01610 canonical_results['model_sum']['value'] = model_sum
01611 canonical_results['model_pbcor_sum'] = {}
01612 canonical_results['model_pbcor_sum']['value'] = model_pbcor_sum
01613 
01614 canonical['results'] = canonical_results
01615 
01616 print "Canonical Regression (default) from "+canonical['date']
01617 
01618 #
01619 # Try and load previous results from regression file
01620 #
01621 regression = {}
01622 regressfile = scriptprefix + '.pickle'
01623 prev_results = {}
01624 
01625 try:
01626     fr = open(regressfile,'r')
01627 except:
01628     print "No previous regression results file "+regressfile
01629     regression['exist'] = False
01630 else:
01631     u = pickle.Unpickler(fr)
01632     regression = u.load()
01633     fr.close()
01634     print "Regression results filled from "+regressfile
01635     print "Regression from version "+regression['version']+" on "+regression['date']
01636     regression['exist'] = True
01637 
01638     prev_results = regression['results']
01639     
01640 #
01641 ##########################################################################
01642 # Calculate regression values
01643 ##########################################################################
01644 #
01645 print '--Calculate Results--'
01646 print ''
01647 #
01648 # Currently using non-PBcor values
01649 #
01650 try:
01651     srcmax = srcstat['max'][0]
01652 except:
01653     srcmax = 0.0
01654 
01655 try:
01656     offrms = offstat['sigma'][0]
01657 except:
01658     offrms = 0.0
01659 
01660 try:
01661     offlinerms = offlinestat['sigma'][0]
01662 except:
01663     offlinerms = 0.0
01664     
01665 try:
01666     momzero_max = momzerostat['max'][0]
01667 except:
01668     momzero_max = 0.0
01669 
01670 try:
01671     momzero_rms = momzerostat['rms'][0]
01672 except:
01673     momzero_rms = 0.0
01674 
01675 try:
01676     momone_median = momonestat['median'][0]
01677 except:
01678     momone_median = 0.0
01679 
01680 #
01681 # Added these sanity checks STM 2008-06-30
01682 #
01683 try:
01684     momone_plane0 = momoneplane0['median'][0]
01685 except:
01686     momone_plane0 = vel0
01687 
01688 try:
01689     momone_plane35 = momoneplane35['median'][0]
01690 except:
01691     momone_plane35 = vel35
01692 
01693 try:
01694     psfcenter = psfstat['mean'][0]
01695 except:
01696     psfcenter = 0.0
01697 
01698 #
01699 # Added model image fluxes STM 2008-10-22
01700 #
01701 try:
01702     modflux = modstat['sum'][0]
01703 except:
01704     modflux = 0.0
01705 
01706 try:
01707     pbcormodflux = pbcormodstat['sum'][0]
01708 except:
01709     pbcormodflux = 0.0
01710 
01711 #
01712 # Store results in dictionary
01713 #
01714 new_regression = {}
01715 
01716 # Some date and version info
01717 import datetime
01718 datestring=datetime.datetime.isoformat(datetime.datetime.today())
01719 
01720 # System info
01721 myvers = casalog.version()
01722 try:
01723     myuser = os.getlogin()
01724 except:
01725     myuser = os.getenv('USER')
01726 #myhost = str( os.getenv('HOST') )
01727 myuname = os.uname()
01728 myhost = myuname[1]
01729 myos = myuname[0]+' '+myuname[2]+' '+myuname[4]
01730 mycwd = os.getcwd()
01731 mypath = os.environ.get('CASAPATH')
01732 
01733 mydataset = 'NGC4826 16apr98 BIMA'
01734 
01735 # Save info in regression dictionary
01736 new_regression['date'] = datestring
01737 new_regression['version'] = myvers
01738 new_regression['user'] = myuser
01739 new_regression['host'] = myhost
01740 new_regression['cwd'] = mycwd
01741 new_regression['os'] = myos
01742 new_regression['uname'] = myuname
01743 new_regression['aipspath'] = mypath
01744 
01745 new_regression['dataset'] = mydataset
01746 
01747 # Fill results
01748 # Note that 'op' tells what to do for the diff :
01749 #    'divf' = abs( new - prev )/prev
01750 #    'diff' = new - prev
01751 
01752 results = {}
01753 
01754 op = 'divf'
01755 tol = 0.08
01756 results['clean_image_max'] = {}
01757 results['clean_image_max']['name'] = 'Clean image max'
01758 results['clean_image_max']['value'] = srcmax
01759 results['clean_image_max']['op'] = op
01760 results['clean_image_max']['tol'] = tol
01761 
01762 results['clean_image_offsrc_max'] = {}
01763 results['clean_image_offsrc_max']['name'] = 'Clean image off-src rms'
01764 results['clean_image_offsrc_max']['value'] = offrms
01765 results['clean_image_offsrc_max']['op'] = op
01766 results['clean_image_offsrc_max']['tol'] = tol
01767 
01768 results['clean_image_offline_max'] = {}
01769 results['clean_image_offline_max']['name'] = 'Clean image off-line rms'
01770 results['clean_image_offline_max']['value'] = offlinerms
01771 results['clean_image_offline_max']['op'] = op
01772 results['clean_image_offline_max']['tol'] = tol
01773 
01774 results['clean_momentzero_max'] = {}
01775 results['clean_momentzero_max']['name'] = 'Moment 0 image max'
01776 results['clean_momentzero_max']['value'] = momzero_max
01777 results['clean_momentzero_max']['op'] = op
01778 results['clean_momentzero_max']['tol'] = tol
01779 
01780 results['clean_momentzero_rms'] = {}
01781 results['clean_momentzero_rms']['name'] = 'Moment 0 image rms'
01782 results['clean_momentzero_rms']['value'] = momzero_rms
01783 results['clean_momentzero_rms']['op'] = op
01784 results['clean_momentzero_rms']['tol'] = tol
01785 
01786 op = 'diff'
01787 tol = 0.1
01788 results['clean_momentone_median'] = {}
01789 results['clean_momentone_median']['name'] = 'Moment 1 image median'
01790 results['clean_momentone_median']['value'] = momone_median
01791 results['clean_momentone_median']['op'] = op
01792 results['clean_momentone_median']['tol'] = 3.0 # km/s.  Was 0.1 before CAS-2163.
01793 
01794 #
01795 # Added these sanity checks STM 2008-06-30
01796 #
01797 results['clean_momentone_planezero'] = {}
01798 results['clean_momentone_planezero']['name'] = 'Moment 1 plane 0'
01799 results['clean_momentone_planezero']['value'] = momone_plane0
01800 results['clean_momentone_planezero']['op'] = op
01801 results['clean_momentone_planezero']['tol'] = tol
01802 
01803 results['clean_momentone_planelast'] = {}
01804 results['clean_momentone_planelast']['name'] = 'Moment 1 plane 35'
01805 results['clean_momentone_planelast']['value'] = momone_plane35
01806 results['clean_momentone_planelast']['op'] = op
01807 results['clean_momentone_planelast']['tol'] = tol
01808 
01809 tol = 0.01
01810 results['clean_psfcenter'] = {}
01811 results['clean_psfcenter']['name'] = 'PSF CH27 at RefPix'
01812 results['clean_psfcenter']['value'] = psfcenter
01813 results['clean_psfcenter']['op'] = op
01814 results['clean_psfcenter']['tol'] = tol
01815 
01816 # And return the ms mean 2008-07-02 STM
01817 op = 'divf'
01818 tol = 0.08
01819 results['vis_mean_cal'] = {}
01820 results['vis_mean_cal']['name'] = 'Vis mean of cal'
01821 results['vis_mean_cal']['value'] = vismean_cal
01822 results['vis_mean_cal']['op'] = op
01823 results['vis_mean_cal']['tol'] = tol
01824 
01825 results['vis_mean_src'] = {}
01826 results['vis_mean_src']['name'] = 'Vis mean of src'
01827 results['vis_mean_src']['value'] = vismean_src
01828 results['vis_mean_src']['op'] = op
01829 results['vis_mean_src']['tol'] = tol
01830 
01831 # model total fluxes 2008-10-22 STM
01832 
01833 results['model_sum'] = {}
01834 results['model_sum']['name'] = 'Model image sum'
01835 results['model_sum']['value'] = modflux
01836 results['model_sum']['op'] = op
01837 results['model_sum']['tol'] = tol
01838 
01839 results['model_pbcor_sum'] = {}
01840 results['model_pbcor_sum']['name'] = 'PBcor Model image sum'
01841 results['model_pbcor_sum']['value'] = pbcormodflux
01842 results['model_pbcor_sum']['op'] = op
01843 results['model_pbcor_sum']['tol'] = tol
01844 
01845 # Now go through and regress
01846 resultlist = ['clean_image_max','clean_image_offsrc_max','clean_image_offline_max',
01847               'clean_momentzero_max','clean_momentzero_rms','clean_momentone_median',
01848               'clean_momentone_planezero','clean_momentone_planelast','clean_psfcenter',
01849               'vis_mean_cal','vis_mean_src','model_sum','model_pbcor_sum']
01850 
01851 for keys in resultlist:
01852     res = results[keys]
01853     prev = None
01854     if prev_results.has_key(keys):
01855         # This is a known regression
01856         prev = prev_results[keys]
01857         results[keys]['test'] = 'Last'
01858     elif canonical_results.has_key(keys):
01859         # Go back to canonical values
01860         prev = canonical_results[keys]
01861         results[keys]['test'] = 'Canon'
01862     if prev:
01863         new_val = res['value']
01864         prev_val = prev['value']
01865         new_diff = new_val - prev_val
01866         if res['op'] == 'divf':
01867             new_diff /= prev_val
01868 
01869         if abs(new_diff) > res['tol']:
01870             new_status = 'Failed'
01871         else:
01872             new_status = 'Passed'
01873         
01874         results[keys]['prev'] = prev_val
01875         results[keys]['diff'] = new_diff
01876         results[keys]['status'] = new_status
01877     else:
01878         # Unknown regression key
01879         results[keys]['prev'] = 0.0
01880         results[keys]['diff'] = 1.0
01881         results[keys]['status'] = 'Missed'
01882         results[keys]['test'] = 'none'
01883 
01884 # Done filling results
01885 new_regression['results'] = results
01886 
01887 # Dataset size info
01888 datasize_raw = 96.0 # MB
01889 datasize_ms = 279.0 # MB
01890 new_regression['datasize'] = {}
01891 new_regression['datasize']['raw'] = datasize_raw
01892 new_regression['datasize']['ms'] = datasize_ms
01893 
01894 #
01895 # Timing
01896 #
01897 if benchmarking:
01898     # Save timing to regression dictionary
01899     new_regression['timing'] = {}
01900 
01901     total = {}
01902     total['wall'] = (endTime - startTime)
01903     total['cpu'] = (endProc - startProc)
01904     total['rate_raw'] = (datasize_raw/(endTime - startTime))
01905     total['rate_ms'] = (datasize_ms/(endTime - startTime))
01906 
01907     new_regression['timing']['total'] = total
01908 
01909     nstages = 15
01910     new_regression['timing']['nstages'] = nstages
01911 
01912     stages = {}
01913     stages[0] = ['import',(import2time-startTime)]
01914     stages[1] = ['concat',(concat2time-import2time)]
01915     stages[2] = ['clearcal',(clearcal2time-concat2time)]
01916     stages[3] = ['listobs',(list2time-clearcal2time)]
01917     stages[4] = ['plotxy',(plotxy2time-list2time)]
01918     stages[5] = ['tflagdata',(flag2time-plotxy2time)]
01919     stages[6] = ['setjy',(setjy2time-flag2time)]
01920     stages[7] = ['gaincal',(gaincal2time-setjy2time)]
01921     stages[8] = ['fluxscale',(fluxscale2time-gaincal2time)]
01922     stages[9] = ['plotcal',(plotcal2time-fluxscale2time)]
01923     stages[10] = ['applycal',(correct2time-plotcal2time)]
01924     stages[11] = ['split',(split2time-correct2time)]
01925     stages[12] = ['clean',(clean2time-split2time)]
01926     stages[13] = ['math/stat',(math2time-clean2time)]
01927     stages[14] = ['moments',(moments2time-math2time)]
01928     
01929     new_regression['timing']['stages'] = stages
01930 
01931 #
01932 ##########################################################################
01933 #
01934 # Save regression results as dictionary using Pickle
01935 #
01936 pickfile = 'out.'+prefix + '.regression.'+datestring+'.pickle'
01937 f = open(pickfile,'w')
01938 p = pickle.Pickler(f)
01939 p.dump(new_regression)
01940 f.close()
01941 
01942 print ""
01943 print "Regression result dictionary saved in "+pickfile
01944 print ""
01945 print "Use Pickle to retrieve these"
01946 print ""
01947 
01948 # e.g.
01949 # f = open(pickfile)
01950 # u = pickle.Unpickler(f)
01951 # clnmodel = u.load()
01952 # polmodel = u.load()
01953 # f.close()
01954 
01955 #
01956 ##########################################################################
01957 #
01958 # Now print out results
01959 # The following writes a logfile for posterity
01960 #
01961 ##########################################################################
01962 #
01963 #outfile='n4826.'+datestring+'.log'
01964 outfile='out.'+prefix+'.'+datestring+'.log'
01965 logfile=open(outfile,'w')
01966 
01967 # Print version info to outfile
01968 print >>logfile,'Running '+myvers+' on host '+myhost
01969 print >>logfile,'at '+datestring
01970 print >>logfile,''
01971 
01972 #
01973 # Report a few key stats
01974 #
01975 print '  NGC4826 Image Cube Max = '+str(srcstat['max'][0])
01976 print "          At ("+str(srcstat['maxpos'][0])+","+str(srcstat['maxpos'][1])+") Channel "+str(srcstat['maxpos'][3])
01977 print '          '+srcstat['maxposf']
01978 print ''
01979 print '          Off-Source Rms = '+str(offstat['sigma'][0])
01980 print '          Signal-to-Noise ratio = '+str(srcstat['max'][0]/offstat['sigma'][0])
01981 print ''
01982 print '          Off-Line   Rms = '+str(offlinestat['sigma'][0])
01983 print '          Signal-to-Noise ratio = '+str(srcstat['max'][0]/offlinestat['sigma'][0])
01984 print ''
01985 
01986 print >>logfile,'  NGC4826 Image Cube Max = '+str(srcstat['max'][0])
01987 print >>logfile,"          At ("+str(srcstat['maxpos'][0])+","+str(srcstat['maxpos'][1])+") Channel "+str(srcstat['maxpos'][3])
01988 print >>logfile,'          '+srcstat['maxposf']
01989 print >>logfile,''
01990 print >>logfile,'          Off-Source Rms = '+str(offstat['sigma'][0])
01991 print >>logfile,'          Signal-to-Noise ratio = '+str(srcstat['max'][0]/offstat['sigma'][0])
01992 print >>logfile,''
01993 print >>logfile,'          Off-Line   Rms = '+str(offlinestat['sigma'][0])
01994 print >>logfile,'          Signal-to-Noise ratio = '+str(srcstat['max'][0]/offlinestat['sigma'][0])
01995 print >>logfile,''
01996 
01997 # Print out comparison:
01998 res = {}
01999 resultlist = ['clean_image_max','clean_image_offsrc_max','clean_image_offline_max',
02000                'clean_momentzero_max','clean_momentzero_rms','clean_momentone_median',
02001               'clean_momentone_planezero','clean_momentone_planelast','clean_psfcenter',
02002               'vis_mean_cal','vis_mean_src','model_sum','model_pbcor_sum']
02003 
02004 # First versus canonical values
02005 print >>logfile,'---'
02006 print >>logfile,'Regression versus previous values:'
02007 print >>logfile,'---'
02008 print '---'
02009 print 'Regression versus previous values:'
02010 print '---'
02011 
02012 if regression['exist']:
02013     print >>logfile,"  Regression results filled from "+regressfile
02014     print >>logfile,"  Regression from version "+regression['version']+" on "+regression['date']
02015     print >>logfile,"  Regression platform "+regression['host']
02016     
02017     print "  Regression results filled from "+regressfile
02018     print "  Regression from version "+regression['version']+" on "+regression['date']
02019     print "  Regression platform "+regression['host']
02020     if regression.has_key('aipspath'):
02021         print >>logfile,"  Regression casapath "+regression['aipspath']
02022         print "  Regression casapath "+regression['aipspath']
02023     
02024 else:
02025     print >>logfile,"  No previous regression file"
02026 
02027 print ""
02028 print >>logfile,""
02029 
02030 final_status = 'Passed'
02031 for keys in resultlist:
02032     res = results[keys]
02033     print '--%30s : %12.6f was %12.6f %4s %12.6f (%6s) %s ' % ( res['name'], res['value'], res['prev'], res['op'], res['diff'], res['status'], res['test'] )
02034     print >>logfile,'--%30s : %12.6f was %12.6f %4s %12.6f (%6s) %s ' % ( res['name'], res['value'], res['prev'], res['op'], res['diff'], res['status'], res['test'] )
02035     if res['status']=='Failed':
02036         final_status = 'Failed'
02037 
02038 if (final_status == 'Passed'):
02039     regstate=True
02040     print >>logfile,'---'
02041     print >>logfile,'Passed Regression test for NGC 4826 Mosaic'
02042     print >>logfile,'---'
02043     print 'Passed Regression test for NGC 4826 Mosaic'
02044     print ''
02045     print 'Regression PASSED'
02046     print ''
02047 else:
02048     regstate=False
02049     print >>logfile,'----FAILED Regression test for NGC 4826 Mosaic'
02050     print '----FAILED Regression test for NGC 4826 Mosaic'
02051     print ''
02052     print 'Regression FAILED'
02053     print ''
02054     
02055 #
02056 ##########################################################################
02057 # Print benchmarking etc.
02058 
02059 if benchmarking:
02060     print ''
02061     print 'Total wall clock time was: %10.3f ' % total['wall']
02062     print 'Total CPU        time was: %10.3f ' % total['cpu']
02063     print 'Raw processing rate MB/s was: %8.1f ' % total['rate_raw']
02064     print 'MS  processing rate MB/s was: %8.1f ' % total['rate_ms']
02065     print ''
02066     print '* Breakdown:                              *'
02067 
02068     print >>logfile,''
02069     print >>logfile,'********* Benchmarking *************************'
02070     print >>logfile,'*                                              *'
02071     print >>logfile,'Total wall clock time was: %10.3f ' % total['wall']
02072     print >>logfile,'Total CPU        time was: %10.3f ' % total['cpu']
02073     print >>logfile,'Raw processing rate MB/s was: %8.1f ' % total['rate_raw']
02074     print >>logfile,'MS  processing rate MB/s was: %8.1f ' % total['rate_ms']
02075     print >>logfile,'* Breakdown:                                   *'
02076 
02077     for i in range(nstages):
02078         print '* %16s * time was: %10.3f ' % tuple(stages[i])
02079         print >>logfile,'* %16s * time was: %10.3f ' % tuple(stages[i])
02080     
02081     print >>logfile,'************************************************'
02082     print >>logfile,'sandrock (2008-06-17) wall time was: 377 seconds'
02083     print >>logfile,'sandrock (2008-06-17) CPU  time was: 312 seconds'
02084 
02085 logfile.close()
02086 
02087 print "Done with NGC4826 Tutorial Regression"
02088 #
02089 ##########################################################################