From 0df7c7b2f184344ee6ecc44c6f9075b168762c25 Mon Sep 17 00:00:00 2001
From: Marco Cammarata <marcocamma@gmail.com>
Date: Thu, 1 Dec 2016 17:02:15 +0100
Subject: [PATCH] more cleanup and add funciton for showing abs spectra with
 vernier sweep

---
 alignment.py             |  1 +
 dispersiveXanes_utils.py | 74 ----------------------------------------
 xppl37_spectra.py        | 49 +++++++++++++++++++-------
 3 files changed, 37 insertions(+), 87 deletions(-)
 create mode 120000 alignment.py

diff --git a/alignment.py b/alignment.py
new file mode 120000
index 0000000..bb969dc
--- /dev/null
+++ b/alignment.py
@@ -0,0 +1 @@
+dispersiveXanes_alignment.py
\ No newline at end of file
diff --git a/dispersiveXanes_utils.py b/dispersiveXanes_utils.py
index e2ace07..10467d7 100644
--- a/dispersiveXanes_utils.py
+++ b/dispersiveXanes_utils.py
@@ -68,77 +68,3 @@ def medianRatio(p1,p2,threshold=0.03):
   return np.ma.average(ratio,axis=0,weights=p1)
 
 
-# /--------------------\
-# |                    |
-# |     PLOTS & CO.    |
-# |                    |
-# \--------------------/
-
-def plotShot(im1,im2,transf1=None,transf2=None,fig=None,ax=None,res=None,E=defaultE,save=None):
-  if transf1 is not None: im1 = transf1.transformImage(im1)
-  if transf2 is not None: im2 = transf2.transformImage(im2)
-  if fig is None and ax is None:
-    fig = plt.subplots(2,3,figsize=[7,5],sharex=True)[0]
-    ax = fig.axes
-  elif fig is not None:
-    ax = fig.axes
-  if E is None: E=np.arange(im1.shape[1])
-  n = im1.shape[0]
-  ax[0].imshow(im1,extent=(E[0],E[-1],0,n),**kw_2dplot)
-  ax[1].imshow(im2,extent=(E[0],E[-1],0,n),**kw_2dplot)
-  ax[2].imshow(im1-im2,extent=(E[0],E[-1],0,n),**kw_2dplot)
-  if res is None:
-    p1 = np.nansum(im1,axis=0)
-    p2 = np.nansum(im2,axis=0)
-    pr = p2/p1
-  else:
-    p1 = res.p1; p2 = res.p2; pr = res.ratio
-  ax[3].plot(E,p1,lw=3)
-  ax[4].plot(E,p1,lw=1)
-  ax[4].plot(E,p2,lw=3)
-  idx = (p1>p1.max()/10.)
-  ax[5].plot(E[idx],pr[idx])
-  if res is not None:
-    ax[5].set_title("FOM: %.2f"%res.fom)
-  else:
-    ax[5].set_title("FOM: %.2f"% calcFOM(p1,p2,pr))
-  if (save is not None) and (save is not False): plt.savefig(save,transparent=True,dpi=500)
-  return fig
-
-def plotRatios(r,shot='random',fig=None,E=defaultE,save=None):
-  if fig is None: fig = plt.subplots(2,1,sharex=True)[0]
-  ax = fig.axes
-  n = r.shape[0]
-  i = ax[0].imshow(r,extent=(E[0],E[-1],0,n),**kw_2dplot)
-  i.set_clim(0,1.2)
-  if shot == 'random' : shot = np.random.random_integers(0,n-1)
-  ax[1].plot(E,r[shot],label="Shot n %d"%shot)
-  ax[1].plot(E,np.nanmedian(r[:10],axis=0),label="median 10 shots")
-  ax[1].plot(E,np.nanmedian(r,axis=0),label="all shots")
-  ax[1].legend()
-  ax[1].set_ylim(0,1.5)
-  ax[1].set_xlabel("Energy")
-  ax[1].set_ylabel("Transmission")
-  ax[0].set_ylabel("Shot num")
-  if (save is not None) and (save is not False): plt.savefig(save,transparent=True,dpi=500)
-
-def plotSingleShots(r,nShots=10,fig=None,E=defaultE,save=None,ErangeForStd=(7090,7150)):
-  if fig is None: fig = plt.subplots(2,1,sharex=True)[0]
-  ax = fig.axes
-  for i in range(nShots):
-    ax[0].plot(E,r[i]+i)
-  ax[0].set_ylim(0,nShots+0.5)
-  av = (1,3,10,30,100)
-  good = np.nanmedian(r,0)
-  for i,a in enumerate(av):
-   m = np.nanmedian(r[:a],0)
-   idx = (E>ErangeForStd[0]) & (E<ErangeForStd[1])
-   fom = np.nanstd( m[idx]/good[idx] ) 
-   print("n shots %d, std %.2f"%(a,fom) )
-   ax[1].plot(E,m+i,label="%d shots, std :%.2f"%(a,fom))
-  ax[1].legend()
-  ax[1].set_ylim(0,len(av)+0.5)
-  ax[1].set_xlabel("Energy")
-  ax[1].set_ylabel("Transmission")
-  if (save is not None) and (save is not False): plt.savefig(save,transparent=True,dpi=500)
-
diff --git a/xppl37_spectra.py b/xppl37_spectra.py
index 0e4dd9e..d0132eb 100644
--- a/xppl37_spectra.py
+++ b/xppl37_spectra.py
@@ -4,6 +4,7 @@ import copy
 import argparse
 import collections
 import dispersiveXanes_alignment as alignment
+import dispersiveXanes_utils as utils
 import xanes_analyzeRun
 
 parser = argparse.ArgumentParser(description='Process argv')
@@ -102,16 +103,16 @@ def calcRef(r1,r2,calibs=None,threshold=0.05):
   out["ratioOfAverage"] = dict()
   out["medianOfRatios"] = dict()
   for p1,p2,n in zip(r1,r2,calibs):
-    out["ratioOfAverage"][n] = alignment.ratioOfAverage(p1,p2,threshold=threshold)
-    out["medianOfRatios"][n] = alignment.medianRatio(p1,p2,threshold=threshold)
+    out["ratioOfAverage"][n] = utils.ratioOfAverage(p1,p2,threshold=threshold)
+    out["medianOfRatios"][n] = utils.medianRatio(p1,p2,threshold=threshold)
   # add curves with all calib together
   p1 = np.vstack(r1)
   p2 = np.vstack(r2)
   n = ",".join(map(str,calibs))
-  ref1 = xanes_analyzeRun.ratioOfAverage(p1,p2,threshold=threshold)
-  ref2 = xanes_analyzeRun.medianRatio(p1,p2,threshold=threshold)
-  out["ratioOfAverage"][n] = alignment.ratioOfAverage(p1,p2,threshold=threshold)
-  out["medianOfRatios"][n] = alignment.medianRatio(p1,p2,threshold=threshold)
+  ref1 = utils.ratioOfAverage(p1,p2,threshold=threshold)
+  ref2 = utils.medianRatio(p1,p2,threshold=threshold)
+  out["ratioOfAverage"][n] = utils.ratioOfAverage(p1,p2,threshold=threshold)
+  out["medianOfRatios"][n] = utils.medianRatio(p1,p2,threshold=threshold)
   out["ratioOfAverage"]['all'] = out["ratioOfAverage"][n]
   out["medianOfRatios"]['all'] = out["medianOfRatios"][n]
   return out
@@ -150,7 +151,7 @@ def calcSampleAbs(run=82,refCalibs=slice(None,None,2),threshold=0.05,refKind="me
   p1 = np.vstack( sample.p1 )
   p2 = np.vstack( sample.p2 )
   print(p1.shape)
-  p1,p2 = alignment.maskLowIntensity(p1,p2,threshold=threshold)
+  p1,p2 = utils.maskLowIntensity(p1,p2,threshold=threshold)
   ratio = p2/p1
   ratio = ratio/ref
   return p1,p2,-np.log10(ratio)
@@ -184,14 +185,27 @@ def showSpectra(run=82,shots=slice(5),calibs=0,averageEachCalib=False,
   ax[0][0].set_xlim(*xlim)
   ax[0][0].set_title("Run %d"%run)
   if not averageEachCalib: ax[0][0].set_ylim(0,shifty*(len(spectra_norm)))
-  
-def showAbs(run=82,shots=slice(5),normalization="auto",shifty=1,xlim=(7080,7180),showAv=True,smoothWidth=0,threshold=0.01):
+
+
+color_ss = '#08519c'
+color_av = '#238b45'
+color_av_all = '#d95f0e'
+
+def showAbs(run=82,shots=slice(5),normalization="auto",shifty=1,xlim=(7080,7180),showAv=True,showAvOverAll=True,smoothWidth=0,threshold=0.01,filterShot=0.1):
   """ normalization: if "auto", the max of the spectra that will be plotted will be used
+      filterShot = means that it filters out the filterShot*100 percentile
   """
   E = alignment.defaultE
   p1,p2,abs = calcSampleAbs(run=run,threshold=threshold)
+  p1_sum = p1.sum(-1)
+  if filterShot>0:
+    idx = p1_sum>np.percentile(p1_sum,filterShot*100)
+    p1 = p1[idx]
+    p2 = p2[idx]
+    abs = abs[idx]
   p1_av  = np.nanmean(p1,axis=0)
   p2_av  = np.nanmean(p2,axis=0)
+  # somehow nanmedian screws up when array is too big ... so using nanmean
   abs_av = np.nanmean(abs,axis=0)
   p1 = p1[shots]; p2=p2[shots]; abs = abs[shots]
   if smoothWidth > 0: abs = smoothSpectra(E,abs,res=smoothWidth)
@@ -204,17 +218,26 @@ def showAbs(run=82,shots=slice(5),normalization="auto",shifty=1,xlim=(7080,7180)
     color = gradual_colors[ishot%len(gradual_colors)]
     ax[0].axhline(ishot*shifty,ls='--',color=color)
     ax[1].axhline(ishot*shifty,ls='--',color=color)
+    if showAvOverAll:
+      if ishot == 0: ax[0].fill_between(E,ishot*shifty,p1_av/normalization+ishot*shifty,color=color_av_all,alpha=0.6)
+      ax[1].plot(E,abs_av+ishot*shifty,color=color_av_all,lw=2,zorder=20)
     if showAv:
-      if ishot == 0: ax[0].fill_between(E,ishot*shifty,p1_av/normalization+ishot*shifty,color='#d95f0e',alpha=0.6)
-      ax[1].plot(E,abs_av+ishot*shifty,color='#d95f0e',lw=2,zorder=20)
+      ax[1].plot(E,np.nanmedian(abs,0)+ishot*shifty,color=color_av,lw=2,zorder=10)
     ax[0].plot(E,s1_norm+ishot*shifty,ls = '-' ,color='0.8',lw=2)
     ax[0].plot(E,s2_norm+ishot*shifty,ls = '-' ,color='0.3',lw=2)
-    ax[1].plot(E,a+ishot*shifty,color=color,lw=2)
+    ax[1].plot(E,a+ishot*shifty,color=color_ss,lw=2)
   ax[0].set_xlim(*xlim)
-  ax[0].set_title("Run %d"%run)
+  ax[0].set_title("Run %s"%str(run))
   ax[1].set_ylabel("Sample Absorption")
   ax[0].set_ylabel("Normalized Spectra")
   ax[0].set_ylim(0,shifty*(p1.shape[0]))
+  print("STD of (average over shown shots) - (average over all): %.3f"%np.nanstd(np.nanmedian(abs,0)-abs_av))
+  ax[1].set_title("STD of (average_{shown}) - (average_{all}): %.3f"%np.nanstd(np.nanmedian(abs,0)-abs_av))
+
+
+def showAbsWithSweep(run=(155,156),first=0,period=150,nSpectra=10,**kwards):
+  shots = slice(first,first+period,int(period/nSpectra))
+  showAbs(run=run,shots=shots,**kwards)
  
 def smoothSpectra(E,abs_spectra,res=0.5):
   from scipy import integrate