fix uneven spectral averaging, update theoretical sensitivity calcula…

…tion and weblog plot

auto_selfcal.py

@@ -1,5 +1,5 @@
 #future improvements
-# get_sensitivity to properly weight the estimated sensitivity by the relative fraction of time on source
+
 # heuristics for switching between calonly and calflag
 # heuristics to switch from combine=spw to combine=''
 # switch heirarchy of selfcal_library such that solint is at a higher level than vis. makes storage of some parameters awkward since they live
@@ -221,8 +221,9 @@
    dirty_SNR,dirty_RMS=estimate_SNR(sani_target+'_'+band+'_dirty.image.tt0')
    dr_mod=1.0
    if telescope =='ALMA' or telescope =='ACA':
-      sensitivity=get_sensitivity(vislist,selfcal_library[target][band][vis]['spws'],spw=selfcal_library[target][band][vis]['spwsarray'],imsize=imsize[band],cellsize=cellsize[band])
+      sensitivity=get_sensitivity(vislist,selfcal_library[target][band],selfcal_library[target][band][vis]['spws'],spw=selfcal_library[target][band][vis]['spwsarray'],imsize=imsize[band],cellsize=cellsize[band])
       dr_mod=get_dr_correction(telescope,dirty_SNR*dirty_RMS,sensitivity,vislist)
+      sensitivity_nomod=sensitivity.copy()
       print('DR modifier: ',dr_mod)
    if not os.path.exists(sani_target+'_'+band+'_initial.image.tt0'):
       if telescope=='ALMA' or telescope =='ACA':
@@ -240,7 +241,7 @@
    initial_NF_SNR,initial_NF_RMS=estimate_near_field_SNR(sani_target+'_'+band+'_initial.image.tt0')
    header=imhead(imagename=sani_target+'_'+band+'_initial.image.tt0')
    if telescope =='ALMA' or telescope == 'ACA':
-      selfcal_library[target][band]['theoretical_sensitivity']=sensitivity
+      selfcal_library[target][band]['theoretical_sensitivity']=sensitivity_nomod
    if 'VLA' in telescope:
       selfcal_library[target][band]['theoretical_sensitivity']=-99.0
    selfcal_library[target][band]['SNR_orig']=initial_SNR
@@ -292,7 +293,7 @@
             dirty_SNR,dirty_RMS=estimate_SNR(sani_target+'_'+band+'_'+spw+'_dirty.image.tt0')
             if not os.path.exists(sani_target+'_'+band+'_'+spw+'_initial.image.tt0'):
                if telescope=='ALMA' or telescope =='ACA':
-                  sensitivity=get_sensitivity(vislist,spw,spw=np.array([int(spw)]),imsize=imsize[band],cellsize=cellsize[band])
+                  sensitivity=get_sensitivity(vislist,selfcal_library[target][band],spw,spw=np.array([int(spw)]),imsize=imsize[band],cellsize=cellsize[band])
                   dr_mod=1.0
                   dr_mod=get_dr_correction(telescope,dirty_SNR*dirty_RMS,sensitivity,vislist)
                   print('DR modifier: ',dr_mod,'SPW: ',spw)
@@ -362,7 +363,7 @@
       selfcal_library[target][band]['nsigma']=10**np.linspace(np.log10(nsigma_init),np.log10(3.0),len(solints[band]))
 
    if telescope=='ALMA' or telescope =='ACA': #or ('VLA' in telescope) 
-      sensitivity=get_sensitivity(vislist,selfcal_library[target][band][vis]['spws'],spw=selfcal_library[target][band][vis]['spwsarray'],imsize=imsize[band],cellsize=cellsize[band])
+      sensitivity=get_sensitivity(vislist,selfcal_library[target][band],selfcal_library[target][band][vis]['spws'],spw=selfcal_library[target][band][vis]['spwsarray'],imsize=imsize[band],cellsize=cellsize[band])
       if band =='Band_9' or band == 'Band_10':   # adjust for DSB noise increase
          sensitivity=sensitivity*4.0 
       if ('VLA' in telescope):
@@ -580,7 +581,7 @@
    vislist=selfcal_library[target][band]['vislist'].copy()
    ## omit DR modifiers here since we should have increased DR significantly
    if telescope=='ALMA' or telescope =='ACA':
-      sensitivity=get_sensitivity(vislist,selfcal_library[target][band][vis]['spws'],spw=selfcal_library[target][band][vis]['spwsarray'],imsize=imsize[band],cellsize=cellsize[band])
+      sensitivity=get_sensitivity(vislist,selfcal_library[target][band],selfcal_library[target][band][vis]['spws'],spw=selfcal_library[target][band][vis]['spwsarray'],imsize=imsize[band],cellsize=cellsize[band])
       dr_mod=1.0
       if not selfcal_library[target][band]['SC_success']: # fetch the DR modifier if selfcal failed on source
          dr_mod=get_dr_correction(telescope,selfcal_library[target][band]['SNR_dirty']*selfcal_library[target][band]['RMS_dirty'],sensitivity,vislist)
@@ -629,7 +630,7 @@
    ## omit DR modifiers here since we should have increased DR significantly
             if not os.path.exists(sani_target+'_'+band+'_'+spw+'_final.image.tt0'):
                if telescope=='ALMA' or telescope =='ACA':
-                  sensitivity=get_sensitivity(vislist,spw,spw=np.array([int(spw)]),imsize=imsize[band],cellsize=cellsize[band])
+                  sensitivity=get_sensitivity(vislist,selfcal_library[target][band],spw,spw=np.array([int(spw)]),imsize=imsize[band],cellsize=cellsize[band])
                   dr_mod=1.0
                   if not selfcal_library[target][band]['SC_success']: # fetch the DR modifier if selfcal failed on source
                      dr_mod=get_dr_correction(telescope,selfcal_library[target][band]['SNR_dirty']*selfcal_library[target][band]['RMS_dirty'],sensitivity,vislist)

selfcal_helpers.py

@@ -760,10 +760,11 @@ def get_SNR_self_update(all_targets,band,vislist,selfcal_library,n_ant,solint_cu
          solint_snr[target][band][solint_next]=selfcal_library[target][band][vislist[0]][solint_curr]['SNR_post']/((n_ant-3)**0.5*(selfcal_library[target][band]['Total_TOS']/solint_float)**0.5)
 
 
-def get_sensitivity(vislist,specmode='mfs',spwstring='',spw=[],chan=0,cellsize='0.025arcsec',imsize=1600,robust=0.5,uvtaper=''):
+def get_sensitivity(vislist,selfcal_library,specmode='mfs',spwstring='',spw=[],chan=0,cellsize='0.025arcsec',imsize=1600,robust=0.5,uvtaper=''):
    sensitivities=np.zeros(len(vislist))
+   TOS=np.zeros(len(vislist))
    counter=0
-   scalefactor=2.5
+   scalefactor=1.0
    for vis in vislist:
       im.open(vis)
       im.selectvis(field='',spw=spwstring)
@@ -791,11 +792,16 @@ def get_sensitivity(vislist,specmode='mfs',spwstring='',spw=[],chan=0,cellsize='
           print('# Data in this spw/MS may be flagged')
           print('#')
           continue
+      #print(sens)
       #print(vis,'Briggs Sensitivity = ', sens[1])
       #print(vis,'Relative to Natural Weighting = ', sens[2])  
       sensitivities[counter]=sens[1]*scalefactor
+      TOS[counter]=selfcal_library[vis]['TOS']
       counter+=1
-   estsens=np.sum(sensitivities)/float(counter)/(float(counter))**0.5
+   #estsens=np.sum(sensitivities)/float(counter)/(float(counter))**0.5
+   #print(estsens)
+   estsens=np.sum(sensitivities*TOS)/np.sum(TOS)
+   print('Estimated Sensitivity: ',estsens)
    return estsens
 
 def LSRKfreq_to_chan(msfile, field, spw, LSRKfreq,spwsarray):
@@ -942,17 +948,23 @@ def flagchannels_from_contdotdat(vis,target,spwsarray):
 
 def get_spw_chanwidths(vis,spwarray):
    widtharray=np.zeros(len(spwarray))
+   bwarray=np.zeros(len(spwarray))
+   nchanarray=np.zeros(len(spwarray))
    for i in range(len(spwarray)):
       tb.open(vis+'/SPECTRAL_WINDOW')
       widtharray[i]=np.abs(np.unique(tb.getcol('CHAN_WIDTH', startrow = spwarray[i], nrow = 1)))
+      bwarray[i]=np.abs(np.unique(tb.getcol('TOTAL_BANDWIDTH', startrow = spwarray[i], nrow = 1)))
+      nchanarray[i]=np.abs(np.unique(tb.getcol('NUM_CHAN', startrow = spwarray[i], nrow = 1)))
       tb.close()
 
-   return widtharray
+   return widtharray,bwarray,nchanarray
 
-def get_spw_chanavg(vis,widtharray,desiredWidth=15.625e6):
+def get_spw_chanavg(vis,widtharray,bwarray,chanarray,desiredWidth=15.625e6):
    avgarray=np.zeros(len(widtharray))
    for i in range(len(widtharray)):
-      avgarray[i]=desiredWidth/widtharray[i]
+      nchan=bwarray[i]/desiredWidth
+      nchan=np.round(nchan)
+      avgarray[i]=chanarray[i]/nchan   
       if avgarray[i] < 1.0:
          avgarray[i]=1.0
    return avgarray
@@ -2178,7 +2190,8 @@ def split_to_selfcal_ms(vislist,band_properties,bands,spectral_average):
           for band in bands:
              desiredWidth=get_desired_width(band_properties[vis][band]['meanfreq'])
              print(band,desiredWidth)
-             band_properties[vis][band]['chan_widths']=get_spw_chanavg(vis,get_spw_chanwidths(vis,band_properties[vis][band]['spwarray']),desiredWidth=desiredWidth)
+             widtharray,bwarray,nchanarray=get_spw_chanwidths(vis,band_properties[vis][band]['spwarray'])
+             band_properties[vis][band]['chan_widths']=get_spw_chanavg(vis,widtharray,bwarray,nchanarray,desiredWidth=desiredWidth)
              print(band_properties[vis][band]['chan_widths'])
              chan_widths=chan_widths+band_properties[vis][band]['chan_widths'].astype('int').tolist()
              if spwstring =='':