From 1f652ab162bc12ad04ebabfbe4a0d82df786da6a Mon Sep 17 00:00:00 2001
From: Marco Cammarata <marcocamma@gmail.com>
Date: Tue, 10 Jan 2017 22:43:22 +0100
Subject: [PATCH] more improvements (including some needed to work after
 background subtraction, see example on salen

---
 xray/azav.py               |   5 +-
 xray/dataReduction.py      |  12 +++--
 xray/example_main_salen.py | 106 +++++++++++++++++++++++++++++++++++++
 xray/id9.py                |  13 +++--
 xray/mask.py               |   5 +-
 xray/storage.py            |   4 +-
 xray/utils.py              |   4 +-
 7 files changed, 134 insertions(+), 15 deletions(-)
 create mode 100644 xray/example_main_salen.py

diff --git a/xray/azav.py b/xray/azav.py
index 4e0cabe..a28039a 100644
--- a/xray/azav.py
+++ b/xray/azav.py
@@ -1,7 +1,8 @@
 from __future__ import print_function,division
 
-import logging as log
-log.basicConfig(level=log.INFO)
+import logging
+log = logging.getLogger(__name__)
+
 
 import numpy as np
 np.seterr(all='ignore')
diff --git a/xray/dataReduction.py b/xray/dataReduction.py
index ab72c28..c960195 100644
--- a/xray/dataReduction.py
+++ b/xray/dataReduction.py
@@ -1,7 +1,8 @@
 from __future__ import print_function,division
 
-import logging as log
-log.basicConfig(level=log.INFO)
+import logging 
+log = logging.getLogger(__name__)
+
 
 import numpy as np
 np.seterr(all='ignore')
@@ -40,7 +41,7 @@ def subtractReferences(i,idx_ref, useRatio = False):
       # normal reference for an on chi, the weighted average
       iref[_i]      = weight_before*ref_before + weight_after*ref_after
     if _i>=idx_ref_after: _ref += 1
-    log.debug("For image %d : %d-%d"%(_i,idx_ref_before,idx_ref_after))
+    log.debug("SubtractRederence For image %d : %d-%d"%(_i,idx_ref_before,idx_ref_after))
   if useRatio:
     i /= iref
   else:
@@ -170,6 +171,9 @@ def errorMask(data,threshold=5):
     # sqrt(len(temp)) = sqrt(numOfDiffs); it is needed to estimate error of single Diff
     idx = np.abs(temp-np.median(temp,axis=0)) > threshold*data.err[iscan]*np.sqrt(len(temp))
     idx_mask.append( idx )
+
+    log.debug("errorMask mask, scanpoint: %s, fraction of q points filtered out (average) %.4e [max %.4e])"%\
+              (data.scan[iscan],idx.sum()/idx.size,max(np.sum(idx,axis=1)/idx.shape[1])) ) 
   if "masks" not in data: data['masks'] = dict()
   data['masks']['error'] = idx_mask
   return data
@@ -188,6 +192,8 @@ def chi2Mask(data,threshold=2):
     # expand along other axis (q ...)
     idx = utils.reshapeToBroadcast(idx,data.diffsInScanPoint[iscan])
     idx_mask.append(idx)
+    log.debug("Chi2 mask, scanpoint: %s, curves filtereout out %d/%d (%.2f%%)"%\
+              (data.scan[iscan],idx.sum(),len(idx),idx.sum()/len(idx)*100) )
   if "masks" not in data: data['masks'] = dict()
   data['masks']['chi2'] = idx_mask
   return data
diff --git a/xray/example_main_salen.py b/xray/example_main_salen.py
new file mode 100644
index 0000000..e130607
--- /dev/null
+++ b/xray/example_main_salen.py
@@ -0,0 +1,106 @@
+from __future__ import print_function,division
+import os
+
+
+# prepare logging
+import logging
+logfname = os.path.splitext(__file__)[0] + ".log"
+logging.basicConfig(level=logging.DEBUG,
+                    format='%(asctime)s %(name)-12s %(levelname)-8s %(message)s',
+                    datefmt='%y-%m-%d %H:%M:%S',
+                    filename=logfname,
+                    filemode='w')
+import numpy as np
+import pylab as plt
+import mcutils as mc
+import mcutils.xray as xray
+from mcutils.xray import id9
+#id9 = xray.id9
+
+from dualtree import dualtree
+
+# use npz files (they can handle more stuff (list of arrays,unicode) than h5py)
+id9.default_extension = '.npz'
+#id9.default_extension = '.h5'
+
+g_default_mask = '../masks/fesalen1_run8_careful.edf'
+
+
+def prepareLog():
+  """ It allows printing to terminal on top of logfile """
+  # define a Handler which writes INFO messages or higher to the sys.stderr
+  console = logging.StreamHandler()
+  console.setLevel(logging.INFO)
+  # set a format which is simpler for console use
+  formatter = logging.Formatter('%(message)s')
+  # tell the handler to use this format
+  console.setFormatter(formatter)
+  # add the handler to the root logger
+  logging.getLogger('').addHandler(console)
+
+prepareLog()
+
+def findCenter():
+  files = xray.utils.getEdfFiles("../fesalen/fesalen1/run8/",nFiles=100)
+  img   = xray.azav.read(files).mean(axis=0) - 10.
+  xray.azav.find_center(img)
+  
+
+def azav(folder,nQ=1500,force=False,saveChi=True,mask=g_default_mask):
+  if isinstance(mask,int):
+    files = xray.utils.getFiles(folder,"*.edf*")
+    img   = xray.azav.pyFAIread(files[0])
+    temp = np.ones_like(img,dtype=bool)
+    temp[:mask] = False
+    mask = temp
+  return id9.doFolder_azav(folder,nQ=nQ,force=force,mask=mask,saveChi=saveChi)
+
+
+def removeBaseline(folder,qlims=(0.6,3),max_iter=30):
+  data = azav(folder)
+  idx = (data.q>qlims[0]) & (data.q<qlims[1])
+  data['q'] = data.q[idx]
+  data['data'] = data.data[:,idx]
+  for i in range(len(data.data)):
+    data['data'][i]=data.data[i]-dualtree.baseline(data.data[i],max_iter=max_iter)
+  fname = os.path.splitext(data.filename)[0] + "_nobaseline" + id9.default_extension
+  data.save(fname)
+  return data
+
+def datared(folder,monitor=(0.5,4),showPlot=True,errMask=5,chi2Mask=2,
+            qlims=(0.6,3),withBaseline=False,storageFile='auto',force=False,**kw):
+  # ice contributes to a lot of noise, filter to q<2
+  if storageFile == 'auto':
+    if withBaseline:
+      storageFile = folder + "/" + "pyfai_1d" + id9.default_extension
+    else:
+      storageFile = folder + "/" + "pyfai_1d_nobaseline" + id9.default_extension
+      if not os.path.exists(storageFile) or force: removeBaseline(folder)
+  data = xray.storage.DataStorage(storageFile)
+  data,diffs = id9.doFolder_dataRed(folder,storageFile=data,monitor=monitor,
+               errMask=errMask,chi2Mask=chi2Mask,qlims=qlims,**kw)
+  if showPlot:
+    xray.utils.plotdiffs(diffs.q,diffs.data,t=diffs.scan,
+         absSignal=diffs.dataAbsAvAll,absSignalScale=30)
+    plt.title(folder + " norm %s" % str(monitor))
+    plt.figure()
+    xray.utils.plotdiffs(diffs.q,diffs.dataAbsAvScanPoint,t=diffs.scan)
+    plt.title(folder + " norm %s" % str(monitor))
+  return data,diffs
+
+
+def doall(folder,force=False,removeBaseline=True):
+  azav(folder,force=force)
+  return datared(folder)
+
+def anaAmplitue(run=6):
+  fname = "../tiox/tiox1/run%d/diffs.npz" % run
+  data  = xray.storage.DataStorage(fname)
+  ranges = ( (1.75,1.85), (2.2,2.4), (3.25,3.4) )
+  nPlot = len(ranges)
+  fig,ax = plt.subplots(nPlot,1,sharex=True)
+  for r,a in zip(ranges,ax):
+    idx = (data.q>r[0]) & (data.q<r[1])
+    amplitude = np.abs(data.data[:,idx]).mean(axis=1)
+    a.plot(data.scan,amplitude,'-o')
+    a.set_title("Range %s"%(str(r)))
diff --git a/xray/id9.py b/xray/id9.py
index a9d915a..cf29abf 100644
--- a/xray/id9.py
+++ b/xray/id9.py
@@ -1,5 +1,5 @@
-import logging as log
-log.basicConfig(level=log.INFO)
+import logging
+log = logging.getLogger(__name__)
 
 import os
 import collections
@@ -47,11 +47,16 @@ def doFolder_azav(folder,nQ=1500,force=False,mask=None,saveChi=True,
 
 def doFolder_dataRed(folder,storageFile='auto',monitor=None,
                      funcForAveraging=np.nanmean,errMask=5,chi2Mask=2,qlims=None):
+  """ storageFile cab the the basename of the file (to look for in folder) or 
+      a DataStorage instance """
 
   if storageFile == 'auto' : storageFile = folder +  "/" + "pyfai_1d" + default_extension
 
-  # read azimuthal averaged curves
-  data = storage.DataStorage(storageFile)
+  if isinstance(storageFile,storage.DataStorage):
+    data = storageFile
+  else:
+    # read azimuthal averaged curves
+    data = storage.DataStorage(storageFile)
 
   if qlims is not None:
     idx = (data.q>qlims[0]) & (data.q<qlims[1])
diff --git a/xray/mask.py b/xray/mask.py
index 1426fa6..5f55e87 100644
--- a/xray/mask.py
+++ b/xray/mask.py
@@ -1,8 +1,9 @@
 from __future__ import print_function
 import sys
 if sys.version_info.major == 2: input=raw_input
-import logging as log
-log.basicConfig(level=log.INFO)
+import logging
+log = logging.getLogger(__name__)
+
 import os
 import numpy as np
 import matplotlib.pyplot as plt
diff --git a/xray/storage.py b/xray/storage.py
index 70c4e70..6ab366b 100644
--- a/xray/storage.py
+++ b/xray/storage.py
@@ -7,8 +7,8 @@ import os
 import h5py 
 import collections
 
-import logging as log
-log.basicConfig(level=log.INFO)
+import logging
+log = logging.getLogger(__name__)  # __name__ is "foo.bar" here
 
 def unwrapArray(a,recursive=True,readH5pyDataset=True):
   """ This function takes an object (like a dictionary) and recursivively
diff --git a/xray/utils.py b/xray/utils.py
index f4e6d44..552c409 100644
--- a/xray/utils.py
+++ b/xray/utils.py
@@ -1,7 +1,7 @@
 from __future__ import print_function,division
 
-import logging as log
-log.basicConfig(level=log.INFO)
+import logging
+log = logging.getLogger(__name__)  # __name__ is "foo.bar" here
 
 import numpy as np
 np.seterr(all='ignore')