Updates to s6 post-processing code

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

postprocess/plot_waveform.py

@@ -170,7 +170,7 @@ for ifo in ifoList:
     plt.xlabel('Time (s)')
     plt.ylabel('Frequency (Hz)')
     plt.axis([axwinner[0], axwinner[1], 0, fs/2])
-    cb = plt.colorbar(cax)
+    cb = plt.colorbar(cax, orientation='horizontal')
     cb.set_label('Logrithmic Signal Power')
     plt.savefig('self_spec_{0}.png'.format(ifo))
     plt.close()

postprocess/s6/many_plots.sh

+#!/bin/bash
+
+# Get path of BayesWaveBurst to figure out where the postprocessing scripts are
+pathtoBWB=$(which BayesWaveBurst)
+
+if [ -s "$pathtoBWB" ] ; then
+    pathtoBWB=${pathtoBWB%src*}
+    fulldirname="$(pwd)"
+    for dirname in $fulldirname/*/*/job*; do
+        echo
+        echo "Entering $dirname ..."
+        cd $dirname
+        MDC="$(grep MDC *bayeswave.run)"
+        echo ${#MDC}
+        if [ ${#MDC} -ne "0" ]
+        then
+            echo
+            echo " * Running plot_waveform..."
+            python -W ignore $pathtoBWB/postprocess/plot_waveform.py --mdc=True
+            #python -W ignore /home/francesco.pannarale/scripts/BayesWaveTools/plot_waveform.py --mdc=True
+            echo
+            echo " * Running plot_evidence..."
+            
+        else
+            echo
+            echo " * Running plot_waveform..."
+            python -W ignore $pathtoBWB/postprocess/plot_waveform.py --mdc=False
+            #python -W ignore /home/francesco.pannarale/scripts/BayesWaveTools/plot_waveform.py --mdc=False
+        fi
+        #echo
+        # echo " * Running skymap..."
+        # $pathtoBWB/postprocess/run_skymap.sh
+        echo
+        echo " * Compressing analyzed data..."
+        tar -cf waveforms.tar waveforms
+        tar -cf chains.tar chains
+        tar -cf snr.tar *snr
+        cd $fulldirname
+        echo
+        echo " * ...DONE!"
+        echo
+    done
+    echo
+    echo
+    echo "To move the postprocessing results in order to make them visible. Run:"
+    echo "  rsync -av --progress <path to this directory> <path to the desired directory in the public_html>  --exclude '*/*snr' --exclude '*/waveforms' --exclude '*/chains' --include='*/' --include '*html' --include '*png' --include condorOut.txt --exclude='*'"
+    echo
+    echo
+else
+    echo "Make sure the BayesWaveBurst command is in your path before running this script!"
+fi

postprocess/s6/plot_evidence.py

-
-
 import matplotlib
 matplotlib.use('Agg')
 import numpy as np
@@ -18,6 +16,10 @@ maxodds=50
 # -- Get directory list
 glitch = []
 signal = []
+scale = []
+print "the type of scale is {0}".format(type(scale))
+print type(scale)
+print type(signal)
 topdir = os.getcwd()
 dirList = glob.glob('*/job*')
 # dirList = [os.path.abspath(directory) for directory in dirList]
@@ -26,12 +28,14 @@ dirList = glob.glob('*/job*')
 # -- Loop over jobs
 snrList = []
 goodJobList = []
+print "the type of scale is {0}".format(type(scale))
 for job in dirList:
     os.chdir(job)
     print "Entered", os.getcwd()
     jobname = glob.glob('*bayeswave.run')[0]
     bayeswave = open(jobname, 'r')
 
+    print "the type of scale is {0}".format(type(scale))
     # ---------------------------------------
     # Extract information from bayeswave.run
     # ---------------------------------------
@@ -39,13 +43,18 @@ for job in dirList:
     for line in bayeswave:
         spl = line.split()
         if len(spl) < 1: continue
-        del spl[-1]
-        # Determine names of IFOs
+        
+        # Determine names of IFOs and MDC scale factor
         if spl[0] == 'Command' and spl[1] == 'line:':
             for index, splElement in enumerate(spl):
                 if splElement == '--ifo':
                     ifoNames.append(spl[index+1]) 
+                if splElement == '--MDC-prefactor':
+                    print spl
+                    scale.append(float(spl[index+1]))
+                    print "At this point, scale is {0}".format(scale)
             continue         
+        del spl[-1]
         # Determine number of IFOs
         if ' '.join(spl) == 'number of detectors':
             spl = line.split()
@@ -74,6 +83,9 @@ for job in dirList:
     else:
         jobname = ''
 
+    # -------------------
+    # Read evidence file
+    # -------------------
     try:
         infile = open(str(jobname)+'evidence.dat', 'r')
     except:
@@ -187,23 +199,21 @@ snrTable.close()
 # ----------------------------
 glitch = np.array(glitch)
 signal = np.array(signal)
+scale = np.array(scale)
+uniq_scale = np.array(list(set(scale)))
 snrList = np.array(snrList)
+eff = np.array([])
+for s in uniq_scale:
+    indx = [scale == s]
+    this_signal = signal[indx]
+    this_glitch = glitch[indx]
+    total_inj = this_signal.size
+    detected = np.logical_and(this_signal>0, this_glitch>0).sum()
+    eff = np.append(eff, detected*1.0/total_inj)
 
-binsize = 2
-snrBins = np.arange(0,20,binsize)
-
-effList = np.array([])
-for bin in snrBins:
-    indx = np.logical_and(bin<snrList, snrList<bin+binsize)
-    this_gl = glitch[indx]
-    this_n  = signal[indx]
-    detected = np.logical_and(this_gl>0, this_n>0)
-    eff = detected.sum()*1.0 / this_gl.size
-    effList = np.append(effList, eff)
-
-plt.figure()
-plt.plot(snrBins + binsize/2.0, effList, 'o')
-plt.xlabel('Network SNR')
+plt.plot(uniq_scale, eff, 'o')
+plt.ylim((0,1))
+plt.xlabel('Scale factor')
 plt.ylabel('Detection Efficiency')
 plt.savefig('eff.png')
 plt.close()
@@ -258,6 +268,8 @@ outfile = open('top.html', 'w')
 outfile.write('<html><body>')
 outfile.write('<h1>BayesWave Burst Injection Run</h1>')
 outfile.write("<h4><a href='./snrTable.html'>Link to table of all injections</a></h4><br/>")
+outfile.write("<h3>Efficiency Curve</h3>")
+outfile.write("<img src='eff.png' alt='eff.png' width=600>")
 outfile.write('<h3>All Injections</h3>')
 outfile.write("<img src='net_snr_all.png' alt='net_snr_all.png' width=600>")
 outfile.write("<h3>All injections, fixed axis scale</h3>")