Prototype bayesogram

git-svn-id: https://svn.ligo.caltech.edu/svn/bayeswave/trunk@136 c56465c9-8126-4a4f-9d7d-ac845eff4865

postprocess/s6/plot_pec.py

@@ -52,6 +52,28 @@ def get_waveform_bigfile(filename, ifo):
             
     return (median_waveform, up_vec, down_vec)
 
+# -- Create the Bayesogram
+def bayesogram(filename, time, timerange=[-0.1,0.1]):
+    strain_master = np.loadtxt(filename)
+    multiply = strain_master.size / time.size
+    X = list(time)*multiply
+    Y = strain_master.flatten()
+    twodrange = [timerange, [Y.min(), Y.max()]]
+    Nx = np.abs( (timerange[1]-timerange[0])/(time[1] - time[0]) )
+    H, xedge, yedge = np.histogram2d(X,Y,bins=[Nx, 100],range=twodrange)
+    return H, xedge, yedge
+    
+
+    for sample in range(median_waveform.size):
+        vector = np.transpose(strain_master)[sample]
+        sort_vec = np.sort(vector)
+        up_quart = sort_vec[ int(0.84*sort_vec.size) ]
+        down_quart = sort_vec[ int(0.16*sort_vec.size) ]
+        down_vec[sample] = down_quart
+        up_vec[sample] = up_quart
+            
+    return (median_waveform, up_vec, down_vec)
+
 
 
 # -- Parse command line arguments
@@ -190,6 +212,30 @@ for ifo in ifoList:
     # -- Save the plot
     plt.savefig('waveform_{0}.png'.format(ifo))
     plt.close()
+
+    # -- Bayesogram 
+    filename = 't_recOut.dat.{0}'.format(ifo)
+    H, xedges, yedges = bayesogram(filename, time, axwinner[0:2])
+
+    # H needs to be rotated and flipped
+    H = np.rot90(H)
+    H = np.flipud(H)
+ 
+    # Mask zeros
+    # Hmasked = np.ma.masked_where(H==0,H) # Mask pixels with a value of zero
+ 
+    # Plot 2D histogram using pcolor
+    fig2 = plt.figure()
+    # plt.pcolormesh(xedges,yedges,Hmasked)
+    plt.pcolormesh(xedges,yedges,H)
+    plt.xlabel('x')
+    plt.ylabel('y')
+    cbar = plt.colorbar()
+    cbar.ax.set_ylabel('Counts')
+    plt.plot(time, median_waveform, 'w', linewidth=1, alpha=0.5)
+    plt.axis(axwinner)
+    plt.savefig('t_bgram_{0}.png'.format(ifo))
+    plt.close()
     
     # -- Attempt to make a spectrogram of median waveform
     ts = time[1] - time[0]