##################################### # # MOPRA SDtasks Regression Script # Position-Switched data # Version STM 2007-03-01 # ##################################### import time import os casapath=os.environ['AIPSPATH'] os.system('rm -rf mopra-* sdregress_mopra* ') #datapath=casapath+'/data/regression/' datapath='/home/sandrock2/smyers/Testing2/Feb07/20070228/mopra-2005-05-08_0350.rpf' copystring='cp -r '+datapath+' .' os.system(copystring) restore() # # NOTE: this tries to reproduce the ASAP tutorial in # http://www.atnf.csiro.au/computing/software/asap/tutorials/tutorial.html # ########################## # # MOPRA Tutorial # Position-Switched data # ########################## startTime=time.time() startProc=time.clock() # Set parameters default('sdcal') infile = 'mopra-2005-05-08_0350.rpf' telescope = 'AT' fluxunit = 'K' specunit = 'channel' calmode = 'quotient' scanaverage = True timeaverage = True polaverage = False tau = 0.0 scanlist = [] iflist = [] kernel = 'none' kwidth = 10 blmode = 'auto' blpoly = 0 edge = [] thresh = 3 avg_limit = 1 masklist = [] sdfile = 'sdregress_mopra.asap' outform = 'asap' plotlevel = 0 sdcal() #Baseline fit: # Scan[0] Beam[0] IF[0] Pol[0] Cycle[0] # p0= -0.223895 # Scan[0] Beam[0] IF[1] Pol[0] Cycle[0] # p0= -0.580879 # # Set rest frequencies # (currently not in task) s = sd.scantable('sdregress_mopra.asap') restfreqs = [110.201,86.243] s.set_restfreqs(restfreqs,"GHz") s.save('sdregress_mopra_fset.asap','ASAP',True) del s # Now do some region statistics # First the line-free region # Set parameters default('sdstat') sdfile = 'sdregress_mopra_fset.asap' telescope = 'AT' fluxunit = 'Jy' specunit = 'km/s' frame = 'LSRK' masklist = [[-70,-30],[30,70]] invertmask = False xstat = {} sdstat() off_stat = xstat print xstat # Then the line region # Set parameters # IF 0 sdfile = 'sdregress_mopra_fset.asap' telescope = 'AT' fluxunit = 'Jy' specunit = 'km/s' frame = 'LSRK' masklist = [[-30,30]] invertmask = False xstat = {} sdstat() line_stat = xstat print xstat # Now do some line fitting # Set parameters # IF 0 default('sdfit') sdfile = 'sdregress_mopra_fset.asap' telescope = 'AT' fluxunit = 'Jy' specunit = 'km/s' frame = 'LSRK' iflist = [0] fitmode = 'list' maskline = [[-50,50]] nfit = [1] fitfile = 'sdregress_mopra_if0.fit' plotlevel = 0 xstat = {} sdfit() fit_if0_stat = xstat print xstat # IF 1 sdfile = 'sdregress_mopra_fset.asap' telescope = 'AT' fluxunit = 'Jy' specunit = 'km/s' frame = 'LSRK' iflist = [1] fitmode = 'list' maskline = [[-50,50]] nfit = [2] fitfile = 'sdregress_mopra_if1.fit' plotlevel = 0 xstat = {} sdfit() fit_if1_stat = xstat print xstat endProc = time.clock() endTime = time.time() exit project = 'Mopra' prolog = 'sdregress.mopra.' # Compare to regression values #-------------------------------------------------- # rms #-------------------------------------------------- #Scan[0] (ORION_SIO) Time[2005/05/08/03:51:25]: # IF[0] = 13.236 #-------------------------------------------------- #Scan[0] (ORION_SIO) Time[2005/05/08/03:51:25]: # IF[1] = 18.642 #-------------------------------------------------- # #-------------------------------------------------- # max #-------------------------------------------------- #Scan[0] (ORION_SIO) Time[2005/05/08/03:51:25]: # IF[0] = 247.960 #-------------------------------------------------- #Scan[0] (ORION_SIO) Time[2005/05/08/03:51:25]: # IF[1] = 326.822 #-------------------------------------------------- # #-------------------------------------------------- # sum #-------------------------------------------------- #Scan[0] (ORION_SIO) Time[2005/05/08/03:51:25]: # IF[0] = 5416.239 #-------------------------------------------------- #Scan[0] (ORION_SIO) Time[2005/05/08/03:51:25]: # IF[1] = 8520.537 #-------------------------------------------------- # #Fitting: #Scan[0] Beam[0] IF[0] Pol[0] Cycle[0] # 0: peak = 234.912 Jy , centre = 7.264 km/s, FWHM = 3.578 km/s # area = 894.756 Jy km/s # #Fitting: #Scan[0] Beam[0] IF[1] Pol[0] Cycle[0] # 0: peak = 234.245 Jy , centre = -5.091 km/s, FWHM = 3.996 km/s # area = 996.349 Jy km/s # 1: peak = 131.814 Jy , centre = 12.587 km/s, FWHM = 5.378 km/s # area = 754.531 Jy km/s # prev_if0_rms=13.236 prev_if0_max=247.960 prev_if0_sum=5416.239 prev_if0_peak_0=234.912 prev_if0_cent_0=7.264 prev_if0_fwhm_0=3.578 prev_if1_rms=18.642 prev_if1_max=326.822 prev_if1_sum=8520.537 prev_if1_peak_0=234.245 prev_if1_cent_0=-5.091 prev_if1_fwhm_0=3.996 prev_if1_peak_1=131.814 prev_if1_cent_1=12.587 prev_if1_fwhm_1=5.378 new_if0_rms=off_stat['rms'][0] new_if0_max=line_stat['max'][0] new_if0_sum=line_stat['sum'][0] if ( fit_if0_stat['nfit'] > 0 ): fitconv_if0=True new_if0_peak_0=fit_if0_stat['peak'][0][0] new_if0_cent_0=fit_if0_stat['cent'][0][0] new_if0_fwhm_0=fit_if0_stat['fwhm'][0][0] err_if0_peak_0=fit_if0_stat['peak'][0][1] err_if0_cent_0=fit_if0_stat['cent'][0][1] err_if0_fwhm_0=fit_if0_stat['fwhm'][0][1] else: fitconv_if0=False new_if0_peak_0=0 new_if0_cent_0=0 new_if0_fwhm_0=0 err_if0_peak_0=0 err_if0_cent_0=0 err_if0_fwhm_0=0 new_if1_rms=off_stat['rms'][1] new_if1_max=line_stat['max'][1] new_if1_sum=line_stat['sum'][1] if ( fit_if1_stat['nfit'] > 0 ): fitconv_if1=True new_if1_peak_0=fit_if1_stat['peak'][0][0] new_if1_cent_0=fit_if1_stat['cent'][0][0] new_if1_fwhm_0=fit_if1_stat['fwhm'][0][0] new_if1_peak_1=fit_if1_stat['peak'][1][0] new_if1_cent_1=fit_if1_stat['cent'][1][0] new_if1_fwhm_1=fit_if1_stat['fwhm'][1][0] err_if1_peak_0=fit_if1_stat['peak'][0][1] err_if1_cent_0=fit_if1_stat['cent'][0][1] err_if1_fwhm_0=fit_if1_stat['fwhm'][0][1] err_if1_peak_1=fit_if1_stat['peak'][1][1] err_if1_cent_1=fit_if1_stat['cent'][1][1] err_if1_fwhm_1=fit_if1_stat['fwhm'][1][1] else: fitconv_if1=False new_if1_peak_0=0 new_if1_cent_0=0 new_if1_fwhm_0=0 new_if1_peak_1=0 new_if1_cent_1=0 new_if1_fwhm_1=0 err_if1_peak_0=0 err_if1_cent_0=0 err_if1_fwhm_0=0 err_if1_peak_1=0 err_if1_cent_1=0 err_if1_fwhm_1=0 diff_if0_rms = abs((prev_if0_rms-new_if0_rms)/prev_if0_rms) diff_if0_max = abs((prev_if0_max-new_if0_max)/prev_if0_max) diff_if0_sum = abs((prev_if0_sum-new_if0_sum)/prev_if0_sum) difpass_if0_rms = (diff_if0_rms<0.05) difpass_if0_max = (diff_if0_max<0.05) difpass_if0_sum = (diff_if0_sum<0.05) difpass_if0 = difpass_if0_rms & difpass_if0_max & difpass_if0_sum if(fitconv_if0): diff_if0_peak_0 = abs((prev_if0_peak_0-new_if0_peak_0)/prev_if0_peak_0) diff_if0_cent_0 = abs((prev_if0_cent_0-new_if0_cent_0)/err_if0_cent_0) diff_if0_fwhm_0 = abs((prev_if0_fwhm_0-new_if0_fwhm_0)/prev_if0_fwhm_0) fitpass_if0_peak_0 = (diff_if0_peak_0<0.05) fitpass_if0_cent_0 = (diff_if0_cent_0<0.25) fitpass_if0_fwhm_0 = (diff_if0_fwhm_0<0.05) fitpass_if0_0 = fitpass_if0_peak_0 & fitpass_if0_cent_0 & fitpass_if0_fwhm_0 fitpass_if0 = fitpass_if0_0 else: diff_if0_peak_0 = 0 diff_if0_cent_0 = 0 diff_if0_fwhm_0 = 0 fitpass_if0 = False diff_if1_rms = abs((prev_if1_rms-new_if1_rms)/prev_if1_rms) diff_if1_max = abs((prev_if1_max-new_if1_max)/prev_if1_max) diff_if1_sum = abs((prev_if1_sum-new_if1_sum)/prev_if1_sum) difpass_if1_rms = (diff_if1_rms<0.05) difpass_if1_max = (diff_if1_max<0.05) difpass_if1_sum = (diff_if1_sum<0.05) difpass_if1 = difpass_if1_rms & difpass_if1_max & difpass_if1_sum if(fitconv_if1): diff_if1_peak_0 = abs((prev_if1_peak_0-new_if1_peak_0)/prev_if1_peak_0) diff_if1_cent_0 = abs((prev_if1_cent_0-new_if1_cent_0)/err_if1_cent_0) diff_if1_fwhm_0 = abs((prev_if1_fwhm_0-new_if1_fwhm_0)/prev_if1_fwhm_0) diff_if1_peak_1 = abs((prev_if1_peak_1-new_if1_peak_1)/prev_if1_peak_1) diff_if1_cent_1 = abs((prev_if1_cent_1-new_if1_cent_1)/err_if1_cent_1) diff_if1_fwhm_1 = abs((prev_if1_fwhm_1-new_if1_fwhm_1)/prev_if1_fwhm_1) fitpass_if1_peak_0 = (diff_if1_peak_0<0.05) fitpass_if1_cent_0 = (diff_if1_cent_0<0.25) fitpass_if1_fwhm_0 = (diff_if1_fwhm_0<0.05) fitpass_if1_peak_1 = (diff_if1_peak_1<0.05) fitpass_if1_cent_1 = (diff_if1_cent_1<0.25) fitpass_if1_fwhm_1 = (diff_if1_fwhm_1<0.05) fitpass_if1_0 = fitpass_if1_peak_0 & fitpass_if1_cent_0 & fitpass_if1_fwhm_0 fitpass_if1_1 = fitpass_if1_peak_1 & fitpass_if1_cent_1 & fitpass_if1_fwhm_1 fitpass_if1 = fitpass_if1_0 & fitpass_if1_1 else: diff_if1_peak_0 = 0 diff_if1_cent_0 = 0 diff_if1_fwhm_0 = 0 diff_if1_peak_1 = 0 diff_if1_cent_1 = 0 diff_if1_fwhm_1 = 0 fitpass_if1 = False if (difpass_if0_max): print '* IF0 Passed spectrum max test ' if (difpass_if0_rms): print '* IF0 Passed spectrum rms test ' if (difpass_if0_sum): print '* IF0 Passed spectrum sum test' if (difpass_if1_max): print '* IF1 Passed spectrum max test ' if (difpass_if1_rms): print '* IF1 Passed spectrum rms test ' if (difpass_if1_sum): print '* IF1 Passed spectrum sum test' if (difpass_if0 & difpass_if1): print '---Passed Stat test for '+project else: print '---FAILED Stat test for '+project print ' ' fitpass=True if (fitconv_if0): if (fitpass_if0_peak_0): print '* IF0 C_0 Passed fit peak test ' if (fitpass_if0_cent_0): print '* IF0 C_0 Passed fit centroid test ' if (fitpass_if0_fwhm_0): print '* IF0 C_0 Passed fit FWHM test' else: print "IF0 failed to converge" fitpass = False if (fitconv_if1): if (fitpass_if1_peak_0): print '* IF1 C_0 Passed fit peak test ' if (fitpass_if1_cent_0): print '* IF1 C_0 Passed fit centroid test ' if (fitpass_if1_fwhm_0): print '* IF1 C_0 Passed fit FWHM test' if (fitpass_if1_peak_1): print '* IF1 C_1 Passed fit peak test ' if (fitpass_if1_cent_1): print '* IF1 C_1 Passed fit centroid test ' if (fitpass_if1_fwhm_1): print '* IF1 C_1 Passed fit FWHM test' else: print "IF0 failed to converge" fitpass = False if (fitpass_if0 & fitpass_if1): print '---Passed Fit test for '+project else: print '---FAILED Fit test for '+project print '' print '' print 'Total wall clock time was: '+str(endTime - startTime) print 'Total CPU time was: '+str(endProc - startProc) #print 'Processing rate MB/s was: ', 35.1/(endTime - startTime) # # NOW TO REGRESSION LOGFILE # import datetime datestring=datetime.datetime.isoformat(datetime.datetime.today()) outfile=prolog+datestring+'.log' logfile=open(outfile,'w') print >>logfile,'' print >>logfile,'************ Regression ****************' print >>logfile,'* *' if (difpass_if0_max): print >>logfile,'* IF0 Passed spectrum max test ' print >>logfile,'* IF0 Spectrum max '+str(new_if0_max) if (difpass_if0_rms): print >>logfile,'* IF0 Passed spectrum rms test ' print >>logfile,'* IF0 Spectrum rms '+str(new_if0_rms) if (difpass_if0_sum): print >>logfile,'* IF0 Passed spectrum sum test' print >>logfile,'* IF0 Spectrum sum '+str(new_if0_sum) if (difpass_if1_max): print >>logfile,'* IF0 Passed spectrum max test ' print >>logfile,'* IF1 Spectrum max '+str(new_if1_max) if (difpass_if1_rms): print >>logfile,'* IF0 Passed spectrum rms test ' print >>logfile,'* IF1 Spectrum rms '+str(new_if1_rms) if (difpass_if1_sum): print >>logfile,'* IF0 Passed spectrum sum test' print >>logfile,'* IF1 Spectrum sum '+str(new_if1_sum) if (difpass_if0 & difpass_if1): print >>logfile,'---' print >>logfile,'Passed Stat test for '+project else: print >>logfile,'---' print >>logfile,'FAILED Stat test for '+project print >>logfile,'* *' print >>logfile,'****************************************' print >>logfile,'* *' if ( fitconv_if0 ): if (fitpass_if0_peak_0): print >>logfile,'* IF0 C_0 Passed fit peak test ' print >>logfile,'* IF0 C_0 Line fit peak '+str(new_if0_peak_0)+' +/- '+str(err_if0_peak_0) if (fitpass_if0_cent_0): print >>logfile,'* IF0 C_0 Passed fit centroid test ' print >>logfile,'* IF0 C_0 Line fit centroid '+str(new_if0_cent_0)+' +/- '+str(err_if0_cent_0) if (fitpass_if0_fwhm_0): print >>logfile,'* IF0 C_0 Passed fit FWHM test' print >>logfile,'* IF0 C_0 Line fit FHWM '+str(new_if0_fwhm_0)+' +/- '+str(err_if0_fwhm_0) else: print >>logfile,'* IF0 Line fit FAILED to converge' if ( fitconv_if1 ): if (fitpass_if1_peak_0): print >>logfile,'* IF1 C_0 Passed fit peak test ' print >>logfile,'* IF1 C_0 Line fit peak '+str(new_if1_peak_0)+' +/- '+str(err_if1_peak_0) if (fitpass_if1_cent_0): print >>logfile,'* IF1 C_0 Passed fit centroid test ' print >>logfile,'* IF1 C_0 Line fit centroid '+str(new_if1_cent_0)+' +/- '+str(err_if1_cent_0) if (fitpass_if1_fwhm_0): print >>logfile,'* IF1 C_0 Passed fit FWHM test' print >>logfile,'* IF1 C_0 Line fit FHWM '+str(new_if1_fwhm_0)+' +/- '+str(err_if1_fwhm_0) if (fitpass_if1_peak_1): print >>logfile,'* IF1 C_1 Passed fit peak test ' print >>logfile,'* IF1 C_1 Line fit peak '+str(new_if1_peak_1)+' +/- '+str(err_if1_peak_1) if (fitpass_if1_cent_1): print >>logfile,'* IF1 C_1 Passed fit centroid test ' print >>logfile,'* IF1 C_1 Line fit centroid '+str(new_if1_cent_1)+' +/- '+str(err_if1_cent_1) if (fitpass_if1_fwhm_1): print >>logfile,'* IF1 C_1 Passed fit FWHM test' print >>logfile,'* IF1 C_1 Line fit FHWM '+str(new_if1_fwhm_1)+' +/- '+str(err_if1_fwhm_1) else: print >>logfile,'* IF1 Line fit FAILED to converge' if ( fitpass_if0 & fitpass_if1 ): print >>logfile,'Passed Fit test for '+project print >>logfile,'---' else: print >>logfile,'---' print >>logfile,'FAILED Fit test for '+project print >>logfile,'****************************************' print >>logfile,'' print >>logfile,'' print >>logfile,'************ Benchmarking **************' print >>logfile,'* *' print >>logfile,'Total wall clock time was: '+str(endTime - startTime) print >>logfile,'Total CPU time was: '+str(endProc - startProc) #print >>logfile,'Processing rate MB/s was: ', 35.1/(endTime - startTime) logfile.close() ########################## # # End Mopra Regression # ##########################