|
casa
$Rev:20696$
|
|
casa
$Rev:20696$
|
00001 ############################################################################# 00002 # # 00003 # Test Name: # 00004 # Regression/Benchmarking Script for NGC 1333 # 00005 # # 00006 # Rationale for Inclusion: # 00007 # This script illustrates data reduction of a mosaic. # 00008 # The dataset is a VLA SiO Mosaic of a YSO in NGC 1333 consisting of # 00009 # 10 source fields in an almost linear array oriented SE-NW. # 00010 # Two spectral windows are centered on red & blue-shifted emission peaks.# 00011 # Central field also observed in 2 spectral widows centered on red and # 00012 # blue-shifted emission peaks. # 00013 # The NGC 1333 Mosaic was observed on 2 separate days: 2may03 & 8may03 # 00014 # The quality of the 2nd dataset is significantly worse and doesn't add # 00015 # much to the image. # 00016 # Observations are of shocked SiO(1-0) emission in a bipolar outflow # 00017 # (~ 43 GHz, VLA Q band receivers). There is no detectable continuum # 00018 # emission in this data. There are 2 overlapping spectral windows with # 00019 # starting frequencies of 43.4197 GHz & 43.4224 GHz. # 00020 # No tipping scans were done during this observation so we cannot do an # 00021 # opacity correction during the gain calibration. # 00022 # Each calibrator is observed in both spectral windows: # 00023 # 0336+323 = Gain calibrator (Sv ~ 2.1 Jy on 2may; 2.8 Jy on 8may) # 00024 # 0542+498 = Flux density calibrator (3c147, Sv = 0.9144 Jy) # 00025 # This flux calibrator is somewhat resolved. you can derive gain # 00026 # solutions with uvrange=[0,50] (details are in the cookbook). # 00027 # 0319+415 = Bandpass calibrator (3c84, Sv ~ 4.9 Jy on 2may; # 00028 # 7.4 Jy on 8may) # 00029 # A good reference antenna for gain & bandpass calibration = ANT 26 for # 00030 # both days. # 00031 # The test data are in uvfits files called N1333_1.UVFITS and # 00032 # N1333_2.UVFITS. # 00033 # # 00034 # The NGC 1333 data on 2 May suffers from some correlator glitches, # 00035 # a correlator re-boot, and some antennas were not fringing. The # 00036 # correlator glitches are easy to find - amplitudes # 00037 # are 100s of Janskys during a glitch. Non-fringing antennas just # 00038 # look like low-amplitude noise and are difficult to find during # 00039 # initial inspection of the data - but gain # 00040 # solution plots show the amplitudes are near 0 for a given # 00041 # timestamp. Even after you get out obvious bad data, you may still # 00042 # want to clip the source data above a threshold. # 00043 # # 00044 # The 8 May data has some non-fringing antennas and is taken under # 00045 # very poor weather conditions. There are major chunks of data on # 00046 # some antennas that were flagged on-line so it is best to delete # 00047 # what remains on these antennas. Clipping of the final calibrated # 00048 # source data is a must. # 00049 # # 00050 # Features Tested: # 00051 # This script illustrates end-to-end processing with CASA as depicted # 00052 # in the following flow-chart. # 00053 # # 00054 # Input Data Process Output Data # 00055 # # 00056 # NGC1333_1.UVFITS---> importuvfits ----> n1333_1.ms # 00057 # | # 00058 # v # 00059 # listobs ----> casapy.log # 00060 # | # 00061 # v # 00062 # flagdata ----> n1333_1.ms # 00063 # | # 00064 # v # 00065 # setjy ----> n1333_1.ms # 00066 # | # 00067 # v # 00068 # gaincal ----> n1333_1.[12].gcal # 00069 # | # 00070 # v # 00071 # bandpass ----> n1333_1.[12].bcal # 00072 # | # 00073 # v # 00074 # fluxscale ----> n1333_1.[12].fcal # 00075 # | # 00076 # v # 00077 # applycal ----> n1333_1.ms # 00078 # | # 00079 # v # 00080 # split ----> n1333_1.src.t2may.ms + # 00081 # | n1333_1.gcal_t[12].2may.ms # 00082 # v # 00083 # NGC1333_2.UVFITS---> importuvfits ----> n1332_2.ms # 00084 # | # 00085 # v # 00086 # listobs ----> casapy.log # 00087 # | # 00088 # v # 00089 # flagdata ----> n1333_2.ms # 00090 # | # 00091 # v # 00092 # setjy ----> n1333_2.ms # 00093 # | # 00094 # v # 00095 # gaincal ----> n1333_2.[12].gcal # 00096 # | # 00097 # v # 00098 # bandpass ----> n1333_2.[12].bcal # 00099 # | # 00100 # v # 00101 # fluxscale ----> n1333_2.[12].fcal # 00102 # | # 00103 # v # 00104 # applycal ----> n1333_2.ms # 00105 # | # 00106 # v # 00107 # split ----> n1333_2.src.t8may.ms + # 00108 # | n1333_2.gcal_t[12].8may.ms # 00109 # v # 00110 # concat ----> n1333_both.ms # 00111 # | # 00112 # v # 00113 # mosaic ----> n1333_both.model + # 00114 # | n1333_both.residual + # 00115 # | n1333_both.image + # 00116 # | n1333_both.flux # 00117 # v # 00118 # immoments ----> n1333_both.src.tom0.red + # 00119 # | n1333_both.src.tom0.blu + # 00120 # | n1333_both.src.tom0.all + # 00121 # | n1333_both.src.tom1.all # 00122 # v # 00123 # Get various statistics # 00124 # # 00125 # Success/failure criteria: # 00126 # Compare statistics with saved results from past trials. Fail # 00127 # when there is a significant variation from previous runs. # 00128 # # 00129 ############################################################################# 00130 # Updated by RRusk 2008-02-08 from tool-based ngc1333_regression.py # 00131 ############################################################################# 00132 # # 00133 # UPDATED BY gmoellen 2008-11-17: # 00134 # # 00135 # There have long been several missteps in the flagging & calibration # 00136 # which I have repaired: # 00137 # # 00138 # o clip of Amp>0.2 was being applied to all fields, and thus spuriously # 00139 # flagging a very large fraction of 0319. # 00140 # o 0319 appears to be highly resolved. Now use 0336 (integrated over # 00141 # the whole observation) instead. (0319 is now dropped entirely) # 00142 # o added use of canned image model for 0542 (3C147_Q.im) # 00143 # o There *is* an opacity tip for 08May! Using opacity=0.062 for 08May # 00144 # and opacity=0.06 for 02May (from tip on 01May) # 00145 # o Use gaincurve=True for *all* calibration solve and apply, not just # 00146 # some. # 00147 # o The phase stability is not too bad, but scan-based solutions # 00148 # (solint='inf') in gaincal are too long, especially for the long # 00149 # scans on the flux density calibrator. Using solint='int' yields a # 00150 # much more consistent set of flux density results for 0336, among the # 00151 # spws, and (marginally) between the two dates. The net effect is to # 00152 # decrease the flux density of 0336 and the ngc1333. # 00153 # o I've added a few a priori flag commands to deal with bad calibrator # 00154 # scans, based on poor solutions and/or calibrated data # 00155 # o The regression tests now do channel-averaging (to avoid goofy edge # 00156 # channels), and the test values have been revised accordingly. The # 00157 # final peak SNR (~9) is better than before (~6), even though the # 00158 # peak flux density is less (due to better f.d. calibration). # 00159 # o Added a variable to optionally disable the imaging/moments portion # 00160 # # 00161 ############################################################################# 00162 import os 00163 import time 00164 import regression_utility as tstutl 00165 00166 # Enable benchmarking? 00167 benchmarking = True 00168 usedasync = False 00169 00170 # Optionally turn off imaging (which is very time-consuming) 00171 doimage = True 00172 00173 # 00174 # Set up some useful variables 00175 # 00176 # Get path to CASA home directory by stipping name from '$CASAPATH' 00177 pathname=os.environ.get('CASAPATH').split()[0] 00178 00179 # This is where the NGC1333 UVFITS data will be 00180 datapath=pathname+'/data/regression/ATST2/NGC1333/' 00181 fitsdata1=datapath+'N1333_1.UVFITS' 00182 fitsdata2=datapath+'N1333_2.UVFITS' 00183 00184 # 3C147 model image for setjy 00185 modelim=pathname+'/data/nrao/VLA/CalModels/3C147_Q.im' 00186 00187 # The testdir where all output files will be kept 00188 testdir='ngc1333_regression' 00189 00190 # The prefix to use for output files. 00191 prefix=testdir+"/"+'n1333_both' 00192 prefix1=testdir+"/"+'n1333_1' 00193 prefix2=testdir+"/"+'n1333_2' 00194 00195 # Make new test directory 00196 # (WARNING! Removes old test directory of the same name if one exists) 00197 tstutl.maketestdir(testdir) 00198 00199 # Start benchmarking 00200 if benchmarking: 00201 startTime = time.time() 00202 startProc = time.clock() 00203 00204 # 00205 #===================================================================== 00206 # 00207 # Import the data from FITS to MS 00208 # 00209 print '*** 02 MAY ***' 00210 print '--Import--' 00211 00212 # Safest to start from task defaults 00213 default('importuvfits') 00214 00215 # Set up the MS filename and save as new global variable 00216 msfile1 = prefix1 + '.ms' 00217 00218 # Use task importuvfits 00219 fitsfile = fitsdata1 00220 vis = msfile1 00221 antnamescheme="new" 00222 importuvfits() 00223 00224 # Record import time 00225 if benchmarking: 00226 importtime = time.time() 00227 00228 # 00229 #===================================================================== 00230 # 00231 # List a summary of the MS 00232 # 00233 print '--Listobs--' 00234 00235 # Don't default this one. Make use of the previous setting of 00236 # vis. Remember, the variables are GLOBAL! 00237 00238 # You may wish to see more detailed information, like the scans. 00239 # In this case use the verbose = True option 00240 verbose = True 00241 00242 listobs() 00243 00244 # Record listing completion time 00245 if benchmarking: 00246 listtime = time.time() 00247 00248 # 00249 #===================================================================== 00250 # Flagging 00251 # 00252 print '--Flagdata--' 00253 00254 # 00255 # The following information on bad data comes from a test report 00256 # created by Debra Shepherd. It is currently available at 00257 # http://aips2.nrao.edu/projectoffice/almatst1.1/TST1.1.data.description.pdf 00258 # 00259 # NGC 1333: bad data on 2 May: 00260 # o Correlator glitches 00261 # - Spwid 1, field 1: 02-May-2003/21:44:00 to 21:50:00 00262 # - Spwid 1, fields 2,3,4,5,6: 02-May-2003/21:40:00 to 21:51:00 00263 # - Spwid 1, field 7: 02-May-2003/21:56:00 to 21:58:00 00264 # - Spwid 2, fields 8,9,10,11,12: 02-May-2003/21:55:00 to 22:20:00 00265 # o Non-fringing antennas: 00266 # - Ant 9 - everything: all spwids, all fields, all times 00267 # - Ant 14 - spwid 1, 02-May-2003/18:36:00 to 18:53:00 00268 # - Ant 14 - spwid 2, 02-May-2003/18:52:00 to 19:10:00 00269 # o Bad bandpass solutions: 00270 # - Ant 22 - spwid 2 (pathological bandpass solution, 00271 # don't trust this antenna). 00272 # o End channels: 00273 # - Flag channels 1,2,3 and 61, 62, 63 (very noisy) 00274 # Flag them to prevent higher noise in some image planes. 00275 # o Even after all this, the calibrated source data still has amplitudes 00276 # that vary from a maximum of 10Jy/channel to 30Jy/channel. 00277 # The uv weighting will properly downweight the poor data. 00278 # 00279 00280 00281 00282 00283 00284 # Use flagdata() in vector mode 00285 default('tflagdata') 00286 vis = msfile1 00287 mode = 'list' 00288 00289 00290 ################################################## 00291 # 00292 # Get rid of the autocorrelations from the MS 00293 # 00294 # Flag antenna with antennaid 8 00295 # 00296 # Flag all data whose amplitude are not in range [0.0,2.0] on the 00297 # parallel hands ( not done ) 00298 # 00299 # Flag data (which is bad) in a time range 00300 # 00301 # Flag all antenna 14, 15 data in the time ranges stated 00302 # 00303 # Sequential flagdata execution: 00304 # 00305 # flagdata(vis=msfile1, autocorr=True, mode='manualflag') 00306 # flagdata(vis=msfile1,antenna='9', mode='manualflag') 00307 # flagdata(vis=msfile1,mode='manualflag',clipexpr='ABS RR', 00308 # clipminmax=[0.0,2.0],clipoutside=True) 00309 # 00310 # flagdata(vis=msfile1,mode='manualflag',clipexpr='ABS LL', 00311 # clipminmax=[0.0,2.0],clipoutside=True) 00312 # flagdata(vis=msfile1,mode='manualflag',timerange='2003/05/02/21:40:58~2003/05/02/22:01:30') 00313 # 00314 # flagdata(vis=msfile1,mode='manualflag', antenna='14', 00315 # timerange='2003/05/02/18:50:50~2003/05/02/19:13:30') 00316 # flagdata(vis=msfile1,mode='manualflag', antenna='15', spw='0', 00317 # timerange='2003/05/02/22:38:49~2003/05/02/22:39:11') 00318 00319 # Parallel flagdata execution: 00320 flagcmds = ["autocorr=True", 00321 "autocorr=False antenna='VA09'", 00322 "autocorr=False timerange='2003/05/02/21:40:58~2003/05/02/22:01:30'", 00323 "autocorr=False antenna='VA14' timerange='2003/05/02/18:50:50~2003/05/02/19:13:30'", 00324 "autocorr=False antenna='VA15' timerange='2003/05/02/22:38:49~2003/05/02/22:39:11' spw='0'"] 00325 00326 inpfile=flagcmds 00327 #autocorr = [true, false , false, false , false ] 00328 #antenna = ['' , 'VA09', '' , 'VA14', 'VA15'] 00329 #timerange = ['' , '' , '2003/05/02/21:40:58~2003/05/02/22:01:30', \ 00330 # '2003/05/02/18:50:50~2003/05/02/19:13:30', \ 00331 # '2003/05/02/22:38:49~2003/05/02/22:39:11'] 00332 #spw = ['' , '' , '' , '' , '0' ] 00333 00334 00335 ################################################### 00336 # Finally, apply all the flag specifications 00337 # 00338 00339 tflagdata() 00340 00341 00342 # Record flagging completion time 00343 if benchmarking: 00344 flagtime = time.time() 00345 00346 00347 # 00348 #===================================================================== 00349 # 00350 # Set the fluxes of the primary calibrator(s) 00351 # 00352 print '--Setjy--' 00353 default('setjy') 00354 00355 setjy(vis=msfile1,field='0542+498_1',modimage=modelim,scalebychan=False,standard='Perley-Taylor 99') 00356 setjy(vis=msfile1,field='0542+498_2',modimage=modelim,scalebychan=False,standard='Perley-Taylor 99') 00357 00358 # Record setjy completion time 00359 if benchmarking: 00360 setjytime = time.time() 00361 00362 # 00363 #===================================================================== 00364 # 00365 # Gain calibration 00366 # 00367 print '--Gaincal--' 00368 default('gaincal') 00369 00370 gtable1 = prefix1 + '.1.gcal' 00371 gaincal(vis=msfile1,caltable=gtable1, 00372 field='0,12,14',spw='0:4~58', gaintype='G', 00373 opacity=0.06,solint='int',combine='',refant='VA27',minsnr=2.,gaincurve=True) 00374 00375 gtable2 = prefix1 + '.2.gcal' 00376 gaincal(vis=msfile1,caltable=gtable2, 00377 field='6,13,15',spw='1:4~58', gaintype='G', 00378 opacity=0.06,solint='int',combine='',refant='VA27',gaincurve=True) 00379 00380 # gaincal calibration completion time 00381 if benchmarking: 00382 gaintime = time.time() 00383 00384 # 00385 #===================================================================== 00386 # 00387 # Bandpass calibration 00388 # 00389 print '--Bandpass--' 00390 default('bandpass') 00391 00392 btable1 = prefix1 + '.1.bcal' 00393 bandpass(vis=msfile1,caltable=btable1, 00394 field='0',spw='0', 00395 opacity=0.06,gaintable=gtable1,interp='nearest', 00396 refant='VA27',solint='inf',combine='scan',gaincurve=True) 00397 btable2 = prefix1 + '.2.bcal' 00398 bandpass(vis=msfile1,caltable=btable2, 00399 field='6',spw='1', 00400 opacity=0.06,gaintable=gtable2,interp='nearest', 00401 refant='VA27',solint='inf',combine='scan',gaincurve=True) 00402 00403 # bandpass calibration completion time 00404 if benchmarking: 00405 bptime = time.time() 00406 00407 # 00408 #===================================================================== 00409 # 00410 # Bootstrap flux scale 00411 # Transfer the flux density from flux calibrater to gain calibraters 00412 # 00413 print '--Fluxscale--' 00414 default('fluxscale') 00415 00416 ftable1 = prefix1 + '.1.fcal' 00417 fluxscale(vis=msfile1,caltable=gtable1,fluxtable=ftable1, 00418 reference='0542+498_1',transfer=['0336+323_1', '0319+415_1']) 00419 ftable2 = prefix1 + '.2.fcal' 00420 fluxscale(vis=msfile1,caltable=gtable2,fluxtable=ftable2, 00421 reference='0542+498_2',transfer=['0336+323_2', '0319+415_2']) 00422 00423 # Record fluxscale completion time 00424 if benchmarking: 00425 fstime = time.time() 00426 00427 # 00428 #===================================================================== 00429 # 00430 # Apply our calibration solutions to the data 00431 # (This will put calibrated data into the CORRECTED_DATA column) 00432 # 00433 print '--ApplyCal--' 00434 default('applycal') 00435 00436 applycal(vis=msfile1, 00437 field='0~5',spw='0', 00438 gaincurve=True,opacity=0.06, 00439 gaintable=[ftable1,btable1], 00440 gainfield='0',calwt=False) 00441 applycal(vis=msfile1, 00442 field='6~11',spw='1', 00443 gaincurve=True,opacity=0.06, 00444 gaintable=[ftable2,btable2], 00445 gainfield='6',calwt=False) 00446 00447 # Record applycal completion time 00448 if benchmarking: 00449 correcttime = time.time() 00450 00451 00452 # 00453 #===================================================================== 00454 # 00455 # Split the target and gain calibrater data 00456 # 00457 print '--Split (target and cals) --' 00458 default('split') 00459 00460 # Split out the corrected data of target 00461 splitms1 = prefix1 + '.src.t2may.ms' 00462 split(vis=msfile1,outputvis=splitms1, 00463 field='1~5,7~11',spw='0;1:0~62', 00464 datacolumn='corrected') 00465 00466 # Record split src data completion time 00467 if benchmarking: 00468 splitsrctime = time.time() 00469 00470 # Split out calibraters 00471 calsplitms1 = prefix1 + '.gcal_t1.2may.ms' 00472 split(vis=msfile1,outputvis=calsplitms1, 00473 field='0',spw='0:0~62', 00474 datacolumn='corrected') 00475 calsplitms2 = prefix1 + '.gcal_t2.2may.ms' 00476 split(vis=msfile1,outputvis=calsplitms2, 00477 field='6',spw='1:0~62', 00478 datacolumn='corrected') 00479 00480 # Record split cal completion time 00481 if benchmarking: 00482 splitcaltime = time.time() 00483 00484 # 00485 #===================================================================== 00486 # 00487 # Import the data from FITS to MS 00488 # 00489 print '*** 08 MAY ***' 00490 print '--Import--' 00491 00492 # Safest to start from task defaults 00493 default('importuvfits') 00494 00495 # Set up the MS filename and save as new global variable 00496 msfile2 = prefix2 + '.ms' 00497 00498 # Use task importuvfits 00499 fitsfile = fitsdata2 00500 vis = msfile2 00501 antnamescheme='new' 00502 importuvfits() 00503 00504 # Note that there will be a ngc5921.ms.flagversions 00505 # containing the initial flags as backup for the main ms flags. 00506 00507 # Record import time 00508 if benchmarking: 00509 importtime2 = time.time() 00510 00511 # 00512 #===================================================================== 00513 # 00514 # List a summary of the MS 00515 # 00516 print '--Listobs--' 00517 00518 # Don't default this one. Make use of the previous setting of 00519 # vis. Remember, the variables are GLOBAL! 00520 00521 # You may wish to see more detailed information, like the scans. 00522 # In this case use the verbose = True option 00523 verbose = False 00524 00525 listobs() 00526 00527 # Record listing completion time 00528 if benchmarking: 00529 listtime2 = time.time() 00530 00531 # 00532 ################################################## 00533 # 00534 # Flagging 00535 # 00536 00537 print '--Flagdata--' 00538 00539 # 00540 # Bad data as identified in report by Debra Shepherd at 00541 # http://aips2.nrao.edu/projectoffice/almatst1.1/TST1.1.data.description.pdf 00542 # 00543 # NGC 1333: bad data on 8 May: 00544 # o Non-fringing antennas: 00545 # - Ant 9 - everything: all spwids, all fields, all times 00546 # - Ant 15 - spwid 1, 08-May-2003/16:32:10 to 16:46:10 00547 # - Ant 14 - spwid 2, 08-May-2003/16:55:10 to 17:10:00 00548 # o Antennas that were mostly flagged on-line - flag the rest of the data: 00549 # - Ants 10 & 19 - everything: all spwids, all fields, all times 00550 # o End channels: 00551 # - Flag channels 1,2,3 and 61, 62, 63 (very noisy) 00552 # o Even after all this, the calibrated source data still has amplitudes that 00553 # vary from a maximum of 40Jy/channel to 140Jy/channel (this data was taken 00554 # under significantly worse conditions than on 2may). I suggest you flag 00555 # the worst timestamps: 00556 # - 08-May-2003/17:03:00 to 17:10:00 00557 # - 08-May-2003/18:50:00 to 18:57:00 00558 # o The uv weighting will properly down-weight the remaining poor data 00559 # (and down-weight the 8may data relative to the better 2may data). 00560 # Now you know why this 2nd day of data doesn't add much to the 00561 # final image RMS. 00562 # 00563 00564 default('tflagdata') 00565 vis = msfile2 00566 mode = 'list' 00567 00568 # 00569 # 00570 ################################################## 00571 # 00572 # Get rid of the autocorrelations from the MS 00573 # 00574 #autocorr[0] = true 00575 00576 # 00577 # Flag all data from antennas 8,9,18 00578 # 00579 00580 #flagdata(vis=msfile2,antenna='9,10,19', mode='manualflag') 00581 00582 #antenna[1] = 'VA09,VA10,VA19' 00583 00584 # 00585 # Flag antenna 14 for timerange. ANTENNAID 14 is ANTENNANAME 15 00586 # 00587 00588 #default('flagdata') 00589 #flagdata(vis=msfile2,mode='manualflag', 00590 # antenna='15', 00591 # timerange='2003/05/08/00:00:00~2003/05/08/20:00:00') 00592 00593 # 00594 # antenna '22' has a bad scan 144 00595 #flagdata(vis=msfile2,mode='manualflag', 00596 # antenna='6,22', 00597 # scan='144') 00598 00599 # 00600 # 00601 # Flag the last data channels 00602 # 00603 00604 #default('flagdata') 00605 #flagdata(vis=msfile2, spw='*:59;60;61;62', mode='manualflag') 00606 00607 #spw[4] = '*:59;60;61;62' 00608 00609 flagcmds = ["autocorr=True", 00610 "antenna='VA09,VA10,VA19'", 00611 "antenna='VA15' timerange='2003/05/08/00:00:00~2003/05/08/20:00:00'", 00612 "antenna='VA06,VA22' scan='144'", 00613 "spw='*:59;60;61;62'"] 00614 00615 inpfile=flagcmds 00616 # 00617 ################################################### 00618 # Finally, apply all the flag specifications 00619 # 00620 00621 tflagdata() 00622 00623 # Record flagging completion time 00624 if benchmarking: 00625 flagtime2 = time.time() 00626 00627 00628 # 00629 #===================================================================== 00630 # 00631 # Set the fluxes of the primary calibrator(s) 00632 # 00633 print '--Setjy--' 00634 default('setjy') 00635 00636 setjy(vis=msfile2,field='0542+498_1',modimage=modelim,scalebychan=False,standard='Perley-Taylor 99') 00637 setjy(vis=msfile2,field='0542+498_2',modimage=modelim,scalebychan=False,standard='Perley-Taylor 99') 00638 00639 # Record setjy completion time 00640 if benchmarking: 00641 setjytime2 = time.time() 00642 00643 # 00644 #===================================================================== 00645 # 00646 # Gain calibration 00647 # 00648 print '--Gaincal--' 00649 default('gaincal') 00650 00651 gtable2_1 = prefix2 + '.1.gcal' 00652 gaincal(vis=msfile2,caltable=gtable2_1, 00653 field='0,12,14',spw='0:4~58', gaintype='G', 00654 opacity=0.062,solint='int',combine='',refant='VA27',gaincurve=True) 00655 gtable2_2 = prefix2 + '.2.gcal' 00656 gaincal(vis=msfile2,caltable=gtable2_2, 00657 field='6,13,15',spw='1:4~58', gaintype='G', 00658 opacity=0.062,solint='int',combine='',refant='VA27',gaincurve=True) 00659 00660 # gaincal calibration completion time 00661 if benchmarking: 00662 gaintime2 = time.time() 00663 00664 00665 # 00666 #===================================================================== 00667 # 00668 # Bandpass calibration 00669 # 00670 print '--Bandpass--' 00671 default('bandpass') 00672 00673 btable2_1 = prefix2 + '.1.bcal' 00674 bandpass(vis=msfile2,caltable=btable2_1, 00675 field='0',spw='0', 00676 opacity=0.062, 00677 gaintable=gtable2_1,interp='nearest', 00678 refant='VA27', 00679 solint='inf',combine='scan',gaincurve=True) 00680 btable2_2 = prefix2 + '.2.bcal' 00681 bandpass(vis=msfile2,caltable=btable2_2, 00682 field='6',spw='1', 00683 opacity=0.062, 00684 gaintable=gtable2_2,interp='nearest', 00685 refant='VA27', 00686 solint='inf',combine='scan',gaincurve=True) 00687 00688 # bandpass calibration completion time 00689 if benchmarking: 00690 bptime2 = time.time() 00691 00692 # 00693 #===================================================================== 00694 # 00695 # Bootstrap flux scale 00696 # Transfer the flux density from flux calibrater to gain calibraters 00697 # 00698 print '--Fluxscale--' 00699 default('fluxscale') 00700 00701 ftable2_1 = prefix2 + '.1.fcal' 00702 fluxscale(vis=msfile2, caltable=gtable2_1, fluxtable=ftable2_1, 00703 reference='0542+498_1',transfer=['0336+323_1', '0319+415_1']) 00704 ftable2_2 = prefix2 + '.2.fcal' 00705 fluxscale(vis=msfile2, caltable=gtable2_2, fluxtable=ftable2_2, 00706 reference='0542+498_2',transfer=['0336+323_2', '0319+415_2']) 00707 00708 # Record fluxscale completion time 00709 if benchmarking: 00710 fstime2 = time.time() 00711 00712 # 00713 #===================================================================== 00714 # 00715 # Apply our calibration solutions to the data 00716 # Do the correction for the above solutions of bandpass and gain, 00717 # including gain curve. 00718 # 00719 print '--ApplyCal--' 00720 default('applycal') 00721 00722 applycal(vis=msfile2, 00723 field='0~5',spw='0', 00724 gaincurve=True,opacity=0.062, 00725 gaintable=[ftable2_1,btable2_1], 00726 gainfield='0',calwt=False) 00727 applycal(vis=msfile2, 00728 field='6~11',spw='1', 00729 gaincurve=True,opacity=0.062, 00730 gaintable=[ftable2_2,btable2_2], 00731 gainfield='6',calwt=False) 00732 00733 00734 # Record applycal completion time 00735 if benchmarking: 00736 correcttime2 = time.time() 00737 00738 # 00739 #===================================================================== 00740 # 00741 # Split the target and gain calibrater data 00742 # 00743 print '--Split (target and cals) --' 00744 default('split') 00745 00746 # Split out the corrected data of target 00747 splitms2 = prefix2 + '.src.t8may.ms' 00748 split(vis=msfile2,outputvis=splitms2, 00749 field='1~5,7~11',spw='0;1:0~62', 00750 datacolumn='corrected') 00751 00752 # Record split src data completion time 00753 if benchmarking: 00754 splitsrctime2 = time.time() 00755 00756 # Split out calibraters 00757 calsplitms2_1 = prefix2 + '.gcal_t1.8may.ms' 00758 split(vis=msfile2,outputvis=calsplitms2_1, 00759 field='0',spw='0:0~62',datacolumn='corrected') 00760 calsplitms2_2 = prefix2 + '.gcal_t2.8may.ms' 00761 split(vis=msfile2,outputvis=calsplitms2_2, 00762 field='6',spw='1:0~62',datacolumn='corrected') 00763 00764 # Record split cal completion time 00765 if benchmarking: 00766 splitcaltime2 = time.time() 00767 00768 00769 # 00770 #===================================================================== 00771 # 00772 # Merge the 2 corrected data sets 00773 # 00774 print '--Concatenate data sets--' 00775 default('concat') 00776 00777 msfileboth = prefix + '.ms' 00778 shellcmd = 'cp -r '+splitms1+' '+msfileboth 00779 os.system(shellcmd) 00780 vis = [msfileboth, splitms2] 00781 concatvis = msfileboth 00782 freqtol='10MHz' 00783 dirtol='1arcsec' 00784 concat() 00785 00786 # Record concatenation completion time 00787 if benchmarking: 00788 concattime = time.time() 00789 00790 # 00791 #===================================================================== 00792 # Imaging using mosaic 00793 # 00794 # 63 channels in each spectral window (spwid) with about 30 channel 00795 # overlap so about 96 total independent channels. Imager will figure 00796 # out how the windows overlap and will place them properly on the grid. 00797 # 00798 # Set total mosaic size with nx, ny; starting with the third channel (2), 00799 # make 18 channels total with 5-channel averaging. The mosaic phase 00800 # centre is fieldid=0. Use both spectral windows. 00801 # 00802 00803 print '--Image data--' 00804 00805 """ 00806 im.open(thems=msfileboth) 00807 im.selectvis(field=range(0,10),spw=[0,1],nchan=63,start=0,step=1) 00808 im.defineimage(nx=800,ny=800,cellx='0.5arcsec',stokes='I',mode='channel',nchan=18,start=2,step=5,phasecenter=0,spw=[0,1]) 00809 ### natural weighting 00810 im.weight(type='natural') 00811 ### Use default primary beam 00812 im.setvp(dovp=True,usedefaultvp=True,dosquint=False) 00813 ### Use the newer mosaic gridder where mosaicing is done in convolution 00814 im.setoptions(padding=1.0,ftmachine='mosaic') 00815 #### Define up to which point in the beam image will be shown 00816 imfluxscale = prefix+'.src.task.fluxscale' 00817 im.setmfcontrol(scaletype='PBCOR',minpb=0.07,constpb=0.4, fluxscale=imfluxscale) 00818 immodel = prefix+'.src.tall.cln.model' 00819 imimage = prefix+'.src.task.image' 00820 imresidual = prefix+'.src.tall.cln.resid' 00821 im.clean(algorithm='mfclark',niter=1,gain=0.1,threshold=5.,model=immodel, 00822 mask='',image=imimage,residual=imresidual,npercycle=30) 00823 im.close() 00824 """ 00825 00826 imagetime=0.0 00827 momenttime=0.0 00828 if doimage: 00829 default('clean') 00830 00831 vis = msfileboth 00832 imagename = prefix 00833 00834 mode = 'channel' 00835 ### image centered on field 0, making 18 channels 00836 ### averaging 5 data channels to 1 image channel 00837 nchan = 18 00838 start = 4 00839 width = 5 00840 psfmode = 'clark' 00841 imsize = [800,800] 00842 cell = ['0.5arcsec','0.5arcsec'] 00843 stokes='I' 00844 #### clean with niter of 1 ...basically doing a just a dirty mosaic 00845 niter = 1 00846 gain = 0.1 00847 threshold = '5.0mJy' 00848 mask = '' 00849 #cleanbox = [] 00850 ###selecting data from 10 fields 0-9 all channels of the 2 spectral windows 00851 field = '0~9' 00852 spw = '0~1' 00853 timerange = '' 00854 restfreq = '' 00855 phasecenter='0' 00856 modelimage = '' 00857 weighting = 'natural' 00858 imagermode='mosaic' 00859 mosweight = False 00860 ftmachine = 'mosaic' 00861 cyclefactor = 1.5 00862 cyclespeedup = -1 00863 scaletype = 'SAULT' 00864 pbcor=False 00865 minpb = 0.1 00866 sigma = '0.001Jy' 00867 targetflux = '1.0Jy' 00868 constrainflux = False 00869 prior = [''] 00870 negcomponent = 2 00871 scales = [0, 3, 10] 00872 interpolation = 'nearest' 00873 async = False 00874 00875 clean() 00876 00877 # Record imager completion time 00878 if benchmarking: 00879 imagetime = time.time() 00880 00881 # 00882 #===================================================================== 00883 # Do moments 00884 # 00885 00886 print '--Calculate moments--' 00887 00888 default('immoments') 00889 00890 imimage=prefix+'.image' 00891 mom0redoutfile = prefix+'.src.tmom0.red' 00892 mom0blueoutfile = prefix+'.src.tmom0.blu' 00893 mom0alloutfile = prefix+'.src.tmom0.all' 00894 mom1alloutfile = prefix+'.src.tmom1.all' 00895 00896 imagename = imimage 00897 moments = [0] 00898 axis = 'spec' 00899 chans = '2~8' 00900 includepix=[0.003,100.0] 00901 excludepix = [-1] 00902 outfile = mom0redoutfile 00903 async = False 00904 immoments() 00905 00906 outfile = mom0blueoutfile 00907 chans='9~15' 00908 immoments() 00909 00910 outfile = mom0alloutfile 00911 chans='2~15' 00912 immoments() 00913 00914 outfile = mom1alloutfile 00915 moments = 1 00916 includepix=[0.02,100.0] 00917 immoments() 00918 00919 # Record moment estimation completion time 00920 if benchmarking: 00921 momenttime = time.time() 00922 00923 """ 00924 ia.open(infile=imimage) 00925 ia.moments(outfile=mom0redoutfile, 00926 moments=0,axis=3,mask='indexin(3,[2:9])',includepix=[0.003,100.0]) 00927 ia.moments(outfile=mom0blueoutfile, 00928 moments=0,axis=3,mask='indexin(3,[9:15])',includepix=[0.003,100.0]) 00929 ia.moments(outfile=mom0alloutfile, 00930 moments=0,axis=3,mask='indexin(3,[2:15])',includepix=[0.003,100.0]) 00931 ia.moments(outfile=mom1alloutfile, 00932 moments=1,axis=3,mask='indexin(3,[2:15])',includepix=[0.02,100.0]) 00933 ia.close() 00934 """ 00935 00936 00937 endProc = time.clock() 00938 endTime = time.time() 00939 00940 # 00941 #===================================================================== 00942 # Do regression test 00943 # 00944 00945 ms.open(thems=splitms1) 00946 ms.selectchannel(1,2,59,1); 00947 src_2may=max(ms.range(items=["amplitude"]).get('amplitude')) 00948 ms.close 00949 00950 ms.open(thems=splitms2) 00951 ms.selectchannel(1,2,59,1); 00952 src_8may=max(ms.range(items=["amplitude"]).get('amplitude')) 00953 ms.close 00954 00955 ms.open(thems=calsplitms2_1) 00956 ms.selectinit(0,T) 00957 ms.selectchannel(1,2,59,1); 00958 gcal1_8may=max(ms.range(items=["amplitude"]).get('amplitude')) 00959 ms.close() 00960 00961 ms.open(thems=calsplitms2_2) 00962 ms.selectchannel(1,2,59,1); 00963 gcal2_8may=max(ms.range(items=["amplitude"]).get('amplitude')) 00964 ms.close() 00965 00966 ms.open(thems=calsplitms1) 00967 ms.selectchannel(1,2,59,1); 00968 gcal1_2may=max(ms.range(items=["amplitude"]).get('amplitude')) 00969 ms.close() 00970 00971 ms.open(thems=calsplitms2) 00972 ms.selectchannel(1,2,59,1); 00973 gcal2_2may=max(ms.range(items=["amplitude"]).get('amplitude')) 00974 ms.close() 00975 00976 thistest_immax=0.0 00977 thistest_imrms=0.0 00978 if doimage: 00979 ia.open(infile=mom0alloutfile) 00980 statistics=ia.statistics() 00981 thistest_immax=statistics['max'][0] 00982 thistest_imrms=statistics['rms'][0] 00983 # 7499.6601 (n1333_both.ms) 00984 00985 # test values 00986 cal1_2may=4.56 # (channel averaged) 00987 cal2_2may=4.00 # (channel averaged) 00988 cal1_8may=5.72 # (channel averaged) 00989 cal2_8may=6.70 # (channel averaged) 00990 src2may=3.13 # (channel averaged) 00991 src8may=8.18 # (channel averaged) 00992 immax=0.694 00993 imrms=0.0729 00994 00995 diff_cal1_2may=abs((cal1_2may-gcal1_2may)/cal1_2may) 00996 diff_cal2_2may=abs((cal2_2may-gcal2_2may)/cal2_2may) 00997 diff_cal1_8may=abs((cal1_8may-gcal1_8may)/cal1_8may) 00998 diff_cal2_8may=abs((cal2_8may-gcal2_8may)/cal2_8may) 00999 diff_src_2may=abs((src2may-src_2may)/src2may) 01000 diff_src_8may=abs((src8may-src_8may)/src8may) 01001 diff_immax=abs((immax-thistest_immax)/immax) 01002 diff_imrms=abs((imrms-thistest_imrms)/imrms) 01003 01004 01005 if not benchmarking: 01006 print '' 01007 print '--- Done ---' 01008 else: 01009 import datetime 01010 datestring=datetime.datetime.isoformat(datetime.datetime.today()) 01011 outfile='ngc1333.'+datestring+'.log' 01012 logfile=open(outfile,'w') 01013 01014 print >>logfile,'' 01015 print >>logfile,'********** Data Summary *********' 01016 print >>logfile,'*********************************' 01017 print >>logfile,'' 01018 print >>logfile,'********** Regression ***********' 01019 print >>logfile,'* *' 01020 regstate = True 01021 if (diff_cal1_2may < 0.05): 01022 print >>logfile,'* Passed cal1 max amplitude test (2may) *' 01023 else: 01024 print >>logfile,'* Failed cal1 max amplitude test (2may) *' 01025 regstate = False 01026 print >>logfile,'* Cal1 max amp (2may) '+str(gcal1_2may)+' ('+str(cal1_2may)+')' 01027 if (diff_cal2_2may < 0.05): 01028 print >>logfile,'* Passed cal2 max amplitude test (2may) *' 01029 else: 01030 print >>logfile,'* Failed cal2 max amplitude test (2may) *' 01031 regstate = False 01032 print >>logfile,'* Cal2 max amp (2may) '+str(gcal2_2may)+' ('+str(cal2_2may)+')' 01033 if (diff_cal1_8may < 0.05): 01034 print >>logfile,'* Passed cal1 max amplitude test (8may) *' 01035 else: 01036 print >>logfile,'* Failed cal1 max amplitude test (8may) *' 01037 regstate = False 01038 print >>logfile,'* Cal1 max amp (8may) '+str(gcal1_8may)+' ('+str(cal1_8may)+')' 01039 if (diff_cal2_8may < 0.05): 01040 print >>logfile,'* Passed cal2 max amplitude test (8may) *' 01041 else: 01042 print >>logfile,'* Failed cal2 max amplitude test (8may) *' 01043 regstate = False 01044 print >>logfile,'* Cal2 max amp (8may) '+str(gcal2_8may)+' ('+str(cal2_8may)+')' 01045 if (diff_src_2may < 0.05): 01046 print >>logfile,'* Passed src max amplitude test (2may) *' 01047 else: 01048 print >>logfile,'* Failed src max amplitude test (2may) *' 01049 regstate = False 01050 print >>logfile,'* Src max amp (2may) '+str(src_2may)+' ('+str(src2may)+')' 01051 if (diff_src_8may < 0.05): 01052 print >>logfile,'* Passed src max amplitude test (8may) *' 01053 else: 01054 print >>logfile,'* Failed src max amplitude test (8may) *' 01055 regstate = False 01056 print >>logfile,'* Src max amp (8may) '+str(src_8may)+' ('+str(src8may)+')' 01057 if (diff_immax < 0.05): 01058 print >>logfile,'* Passed image max test *' 01059 else: 01060 if doimage: 01061 print >>logfile,'* Failed image max test *' 01062 regstate = False 01063 else: 01064 print >>logfile,'* Did not do image max test' 01065 01066 print >>logfile,'* Image max '+str(thistest_immax)+' ('+str(immax)+')' 01067 if (diff_imrms < 0.05): 01068 print >>logfile,'* Passed image rms test *' 01069 else: 01070 if doimage: 01071 print >>logfile,'* Failed image rms test *' 01072 regstate = False 01073 else: 01074 print >>logfile,'* Did not do image rms test' 01075 print >>logfile,'* Image rms '+str(thistest_imrms)+' ('+str(imrms)+')' 01076 #if ((diff_cal1_2may<0.05) & (diff_cal2_2may<0.05) & 01077 # (diff_cal1_8may<0.05) & (diff_cal2_8may<0.05) & 01078 # (diff_src_2may<0.05) & (diff_src_8may<0.05) & 01079 # (diff_immax<0.05) & (diff_imrms<0.05)): 01080 if regstate: 01081 print >>logfile,'---' 01082 print >>logfile,'Passed Regression test for NGC1333' 01083 print >>logfile,'---' 01084 tstutl.note("Passed Regression test for NGC1333","NORMAL") 01085 else: 01086 print >>logfile,'----FAILED Regression test for NGC1333' 01087 tstutl.note("FAILED Regression test for NGC1333","SEVERE") 01088 print >>logfile,'*********************************' 01089 01090 print >>logfile,'' 01091 print >>logfile,'' 01092 print >>logfile,'********* Benchmarking *****************' 01093 print >>logfile,'* *' 01094 print >>logfile,'Total wall clock time was: '+str(endTime - startTime) 01095 print >>logfile,'Total CPU time was: '+str(endProc - startProc) 01096 print >>logfile,'Processing rate MB/s was: '+str(240.3/(endTime - startTime)) 01097 print >>logfile,'* Breakdown: *' 01098 print >>logfile,'* import time was: '+str(importtime-startTime) 01099 print >>logfile,'* listobs time was: '+str(listtime-importtime) 01100 print >>logfile,'* flagdata time was: '+str(flagtime-listtime) 01101 print >>logfile,'* setjy time was: '+str(setjytime-flagtime) 01102 print >>logfile,'* gaincal time was: '+str(gaintime-setjytime) 01103 print >>logfile,'* bandpass time was: '+str(bptime-gaintime) 01104 print >>logfile,'* fluxscale time was: '+str(fstime-bptime) 01105 print >>logfile,'* correct time was: '+str(correcttime-fstime) 01106 print >>logfile,'* split time was: '+str(splitcaltime-correcttime) 01107 print >>logfile,'* import time was: '+str(importtime2-splitcaltime) 01108 print >>logfile,'+ listobs time was: '+str(listtime2-importtime2) 01109 print >>logfile,'* flagdata time was: '+str(flagtime2-listtime2) 01110 print >>logfile,'* setjy time was: '+str(setjytime2-flagtime2) 01111 print >>logfile,'* gaincal time was: '+str(gaintime2-setjytime2) 01112 print >>logfile,'* bandpass time was: '+str(bptime2-gaintime2) 01113 print >>logfile,'* fluxscale time was: '+str(fstime2-bptime2) 01114 print >>logfile,'* correct time was: '+str(correcttime2-fstime2) 01115 print >>logfile,'* split time was: '+str(splitcaltime2-correcttime2) 01116 print >>logfile,'* concatenate time was: '+str(concattime-splitcaltime2) 01117 print >>logfile,'* image time was: '+str(imagetime-concattime) 01118 print >>logfile,'* moments time was: '+str(momenttime-imagetime) 01119 print >>logfile,'*****************************************' 01120 # 01121 logfile.close()
1.8.0
