new submodules (peaks,cell) and new functions (like backgorundSubtraction (in azav.py)

xray/__init__.py

@@ -1,6 +1,7 @@
 from . import storage
 from . import utils
 from . import mask
+from . import peaks
 try:
   from . import azav
 except ImportError:

xray/azav.py

@@ -52,7 +52,14 @@ def ai_as_dict(ai):
 
 def ai_as_str(ai):
   """ ai is a pyFAI azimuthal intagrator"""
-  return "#" + str(ai).replace("\n","\n#")
+  s=[ "# Detector        : %s" % ai.detector.name,
+      "# Pixel      [um] : %.2fx%.2f" % (ai.pixel1*1e6,ai.pixel2*1e6),
+      "# Distance   [mm] : %.3f" % (ai.dist*1e3),
+      "# Center     [mm] : %.3f,%.3f" % (ai.poni1*1e3,ai.poni2*1e3),
+      "# Center     [px] : %.3f,%.3f" % (ai.poni1/ai.pixel1,ai.poni2/ai.pixel2),
+      "# Wavelength [A]  : %.5f" % (ai.wavelength*1e10),
+      "# rot[1,2,3] [rad]: %.3f,%.3f,%.3f" % (ai.rot1,ai.rot2,ai.rot3) ]
+  return "\n".join(s)
 
 def do1d(ai, imgs, mask = None, npt_radial = 600, method = 'csr',safe=True,dark=10., polCorr = 1):
     """ ai is a pyFAI azimuthal intagrator 
@@ -218,7 +225,9 @@ def leastsq_circle(x,y):
   residu   = np.sum((Ri - R)**2)
   return xc, yc, R
 
-def find_center(img,psize=100e-6,dist=0.1,wavelength=0.8e-10,**kwargs):
+def find_center(img,psize=100e-6,dist=0.1,wavelength=0.8e-10,center=None,reference=None,**kwargs):
+  """ center is the initial centr (can be None)
+      reference is a reference position to be plot in 2D plots """
   plt.ion()
   kw = dict( pixel1 = psize, pixel2 = psize, dist = dist,wavelength=wavelength )
   kw.update(kwargs)
@@ -232,23 +241,30 @@ def find_center(img,psize=100e-6,dist=0.1,wavelength=0.8e-10,**kwargs):
   ans = ""
   print("Enter 'end' when done")
   while ans != "end":
-    if ans == "":
+    if center is None:
       print("Click on beam center:")
       plt.sca(ax_img); # set figure to use for mouse interaction
-      xc,yc = plt.ginput()[0]
+      center = plt.ginput()[0]
-    else:
+    print("Selected center:",center)
-      xc,yc = map(float,ans.split(","))
+    ai.set_poni1(center[0]*psize)
-    print("Selected center:",xc,yc)
+    ai.set_poni2(center[1]*psize)
-    ai.set_poni1(xc*psize)
-    ai.set_poni2(yc*psize)
     q,az,i = do2d(ai,img)
     mesh = ax_pyfai.pcolormesh(q,az,i)
     mesh.set_clim( *np.percentile(i,(2,95) ) )
-    ax_pyfai.set_title(str( (xc,yc) ))
+    ax_pyfai.set_title(str(center))
+    if reference is not None: ax_pyfai.axvline(reference)
     plt.pause(0.01)
     plt.draw()
     ans=input("Enter to continue with clinking or enter xc,yc values ")
-  print("Final values: (in pixels) %.3f %.3f"%(xc,yc))
+    if ans == '':
+      center = None
+    else:
+      try:
+        center = list(map(float,ans.split(",")))
+      except Exception as e:
+        center = None
+    if center == []: center = None
+  print("Final values: (in pixels) %.3f %.3f"%(center[0],center[1]))
   return ai
 
 def average(fileOrFolder,delays=slice(None),scale=1,norm=None,returnAll=False,plot=False,

xray/cell.py

@@ -0,0 +1,68 @@
+from __future__ import division,print_function
+import numpy as np
+from numpy import sin,cos
+
+class Triclinic(object):
+  def __init__(self,a=1,b=1,c=1,alpha=90,beta=90,gamma=90):
+    self.a = a
+    self.b = b
+    self.c = c
+    alpha  = alpha*np.pi/180
+    beta   = beta*np.pi/180
+    gamma  = gamma*np.pi/180
+    self.alpha = alpha
+    self.beta  = beta 
+    self.gamma = gamma
+
+    self._s11 = b**2 * c**2 * sin(alpha)**2
+    self._s22 = a**2 * c**2 * sin(beta)**2
+    self._s33 = a**2 * b**2 * sin(gamma)**2
+    self._s12 = a*b*c**2*(cos(alpha) * cos(beta) - cos(gamma))
+    self._s23 = a**2*b*c*(cos(beta) * cos(gamma) - cos(alpha))
+    self._s13 = a*b**2*c*(cos(gamma) * cos(alpha) - cos(beta))
+    self.V    = (a*b*c)*np.sqrt(1-cos(alpha)**2 - cos(beta)**2 - cos(gamma)**2 + 2*cos(alpha)*cos(beta)*cos(gamma))
+
+  def __call__(self,h,k,l):  return self.q(h,k,l)
+
+  def d(self,h,k,l):
+    temp = self._s11*h**2 + \
+           self._s22*k**2 + \
+           self._s33*l**2 + \
+           2*self._s12*h*k+ \
+           2*self._s23*k*l+ \
+           2*self._s13*h*l
+    d    = self.V/np.sqrt(temp)
+    return d
+
+  def q(self,h,k,l):
+    return 2*np.pi/self.d(h,k,l)
+
+class Orthorombic(Triclinic):
+  def __init__(self,a=1,b=1,c=1):
+    Triclinic.__init__(self,a=a,b=b,c=c,alpha=90,beta=90,gamma=90)
+
+class Monoclinic(object):
+  def __init__(self,a=1,b=1,c=1,beta=90.):
+    self.a = a
+    self.b = b
+    self.c = c
+    beta   = beta/np.pi*180
+    self.beta = beta
+
+    self.V    = (a*b*c)
+
+  def __call__(self,h,k,l):  return self.Q(h,k,l)
+
+  def Q(self,h,k,l):
+    temp = h**2/self.a**2 + (k*sin(self.beta))**2/self.b**2+l**2/self.c**2+2*h*l*cos(self.beta)/self.a/self.c
+    d    = 1/np.sqrt(temp)
+    print(d)
+    return 2*np.pi/d
+
+
+
+ti3o5_lambda = Triclinic(a = 9.83776, b = 3.78674, c = 9.97069, beta = 91.2567)
+ti3o5_beta   = Triclinic(a = 9.7382 , b = 3.8005 , c = 9.4333 , beta = 91.496)
+#ti3o5_beta   = Monoclinic(a = 9.7382 , b = 3.8005 , c = 9.4333 , beta = 91.496)
+ti3o5_alpha  = Triclinic(a = 9.8372,  b = 3.7921,  c = 9.9717)
+#ti3o5_alpha1  = Orthorombic(a = 9.8372,  b = 3.7921,  c = 9.9717)

xray/dataReduction.py

@@ -73,8 +73,10 @@ def averageScanPoints(scan,data,isRef=None,lpower=None,useRatio=False,\
   if isRef.sum()>0:
     # create a copy (subtractReferences works in place)
     diff = subtractReferences(data.copy(),np.argwhere(isRef), useRatio=useRatio)
+    avNeg = funcForAveraging(data[isRef],axis=0)
   else:
     diff = data
+    avNeg = np.zeros_like(avData)
 
   # normalize signal for laser intensity if provided
   if lpower is not None:
@@ -116,10 +118,10 @@ def averageScanPoints(scan,data,isRef=None,lpower=None,useRatio=False,\
 
     # store error of mean
     err[i] = noise/np.sqrt(shot_idx.sum())
-
+  ret = dict(scan=scan_pos,data=ret,dataUnmasked=ret.copy(),err=err,
-  ret = dict(scan=scan_pos,data=ret,dataUnmasked=ret.copy(),err=err,errUnmasked=err.copy(),
+        errUnmasked=err.copy(),chi2_0=chi2_0,diffsInScanPoint=diffsInScanPoint,
-        chi2_0=chi2_0,diffsInScanPoint=diffsInScanPoint,dataAbsAvAll=avData,dataAbsAvScanPoint=data_abs,
+        dataAbsAvNeg = avNeg, dataAsAbs=ret+avNeg,
-        dataAbs=data.copy())
+        dataAbsAvAll=avData,dataAbsAvScanPoint=data_abs,dataAbs=data.copy())
   ret = storage.DataStorage(ret)
   return ret
 

xray/dspace.py

@@ -0,0 +1,70 @@
+from __future__ import division,print_function
+import numpy as np
+from numpy import sin,cos
+
+
+
+class Triclinic(object):
+  def __init__(self,a=1,b=1,c=1,alpha=90,beta=90,gamma=90):
+    self.a = a
+    self.b = b
+    self.c = c
+    alpha  = alpha*np.pi/180
+    beta   = beta*np.pi/180
+    gamma  = gamma*np.pi/180
+    self.alpha = alpha
+    self.beta  = beta 
+    self.gamma = gamma
+
+    self._s11 = b**2 * c**2 * sin(alpha)**2
+    self._s22 = a**2 * c**2 * sin(beta)**2
+    self._s33 = a**2 * b**2 * sin(gamma)**2
+    self._s12 = a*b*c**2*(cos(alpha) * cos(beta) - cos(gamma))
+    self._s23 = a**2*b*c*(cos(beta) * cos(gamma) - cos(alpha))
+    self._s13 = a*b**2*c*(cos(gamma) * cos(alpha) - cos(beta))
+    self.V    = (a*b*c)*np.sqrt(1-cos(alpha)**2 - cos(beta)**2 - cos(gamma)**2 + 2*cos(alpha)*cos(beta)*cos(gamma))
+
+  def __call__(self,h,k,l):  return self.q(h,k,l)
+
+  def d(self,h,k,l):
+    temp = self._s11*h**2 + \
+           self._s22*k**2 + \
+           self._s33*l**2 + \
+           2*self._s12*h*k+ \
+           2*self._s23*k*l+ \
+           2*self._s13*h*l
+    d    = self.V/np.sqrt(temp)
+    return d
+
+  def q(self,h,k,l):
+    return 2*np.pi/self.d(h,k,l)
+
+class Orthorombic(Triclinic):
+  def __init__(self,a=1,b=1,c=1):
+    Triclinic.__init__(self,a=a,b=b,c=c,alpha=90,beta=90,gamma=90)
+
+class Monoclinic(object):
+  def __init__(self,a=1,b=1,c=1,beta=90.):
+    self.a = a
+    self.b = b
+    self.c = c
+    beta   = beta/np.pi*180
+    self.beta = beta
+
+    self.V    = (a*b*c)
+
+  def __call__(self,h,k,l):  return self.Q(h,k,l)
+
+  def Q(self,h,k,l):
+    temp = h**2/self.a**2 + (k*sin(self.beta))**2/self.b**2+l**2/self.c**2+2*h*l*cos(self.beta)/self.a/self.c
+    d    = 1/np.sqrt(temp)
+    print(d)
+    return 2*np.pi/d
+
+
+
+ti3o5_lambda = Triclinic(a = 9.83776, b = 3.78674, c = 9.97069, beta = 91.2567)
+ti3o5_beta   = Triclinic(a = 9.7382 , b = 3.8005 , c = 9.4333 , beta = 91.496)
+#ti3o5_beta   = Monoclinic(a = 9.7382 , b = 3.8005 , c = 9.4333 , beta = 91.496)
+ti3o5_alpha  = Triclinic(a = 9.8372,  b = 3.7921,  c = 9.9717)
+#ti3o5_alpha1  = Orthorombic(a = 9.8372,  b = 3.7921,  c = 9.9717)

xray/id9.py

@@ -45,18 +45,21 @@ def doFolder_azav(folder,nQ=1500,force=False,mask=None,saveChi=True,
   return azav.doFolder(folder,files="*.edf*",nQ=nQ,force=force,mask=mask,
     saveChi=saveChi,poni=poni,storageFile=storageFile,diagnostic=diag)
 
-def doFolder_dataRed(folder,storageFile='auto',monitor=None,
+def doFolder_dataRed(azavStorage,monitor=None,funcForAveraging=np.nanmean,
-                     funcForAveraging=np.nanmean,errMask=5,chi2Mask=2,qlims=None):
+                     errMask=5,chi2Mask=2,qlims=None,outStorageFile='auto'):
-  """ storageFile cab the the basename of the file (to look for in folder) or 
+  """ azavStorage if a DataStorage instance or the filename to read """
-      a DataStorage instance """
 
-  if storageFile == 'auto' : storageFile = folder +  "/" + "pyfai_1d" + default_extension
+  if isinstance(azavStorage,storage.DataStorage):
-
+    data = azavStorage
-  if isinstance(storageFile,storage.DataStorage):
+    folder = os.path.dirname(data.filename) if data.filename is not None else "./"
-    data = storageFile
+  elif os.path.exists(azavStorage):
+    folder = os.path.dirname(azavStorage)
+    data = storage.DataStorage(azavStorage)
   else:
-    # read azimuthal averaged curves
+    # assume is just a folder name
-    data = storage.DataStorage(storageFile)
+    folder = azavStorage
+    azavStorage  = folder +  "/pyfai_1d" + default_extension
+    data = storage.DataStorage(azavStorage)
 
   if qlims is not None:
     idx = (data.q>qlims[0]) & (data.q<qlims[1])
@@ -76,6 +79,8 @@ def doFolder_dataRed(folder,storageFile='auto',monitor=None,
   
   # save txt and npz file
   dataReduction.saveTxt(folder,diffs,info=data.pyfai_info)
-  diffs.save(folder + "/" + "diffs" + default_extension)
+  if outStorageFile == 'auto':
+    outStorageFile = folder + "/diffs" + default_extension
+  diffs.save(outStorageFile)
 
   return data,diffs

xray/peaks.py

@@ -0,0 +1,20 @@
+import lmfit
+import numpy as np
+
+pv = lmfit.models.PseudoVoigtModel()
+
+def fitPeak(x,y,autorange=False):
+  if autorange:
+    # find fwhm
+    idx = np.ravel(np.argwhere( y<y.max()/2 ))
+    # find first crossing
+    p1 = idx[idx<np.argmax(y)][-1]
+    p2 = idx[idx>np.argmax(y)][0]
+    c = int( (p1+p2)/2 )
+    dp = int( np.abs(p1-p2) )
+    idx = slice(c-dp,c+dp)
+    x = x[idx]
+    y = y[idx]
+  pars = pv.guess(y,x=x)
+  ret  = pv.fit(y,x=x,params=pars)
+  return ret

xray/storage.py

@@ -27,7 +27,7 @@ def unwrapArray(a,recursive=True,readH5pyDataset=True):
       a = dict(a); # convert to dict, otherwise can't asssign values
       for key,value in a.items(): a[key] = unwrapArray(value)
     elif isinstance(a,list):
-      for index in range(len(a)): a[index] = unwrapArray(a[i])
+      for index in range(len(a)): a[index] = unwrapArray(a[index])
     else:
       pass
   if isinstance(a,dict): a = DataStorage(a)
@@ -80,25 +80,29 @@ def read(fname):
   extension = os.path.splitext(fname)[1]
   log.info("Reading storage file %s"%fname)
   if extension == ".npz": 
-    return npzToDict(fname)
+    return DataStorage(npzToDict(fname))
   elif extension == ".h5":
-    return h5ToDict(fname)
+    return DataStorage(h5ToDict(fname))
   else:
     raise ValueError("Extension must be h5 or npz, it was %s"%extension) 
 
 def save(fname,d):
   extension = os.path.splitext(fname)[1]
   log.info("Saving storage file %s"%fname)
-  if extension == ".npz": 
+  try:
-    return dictToNpz(fname,d)
+    if extension == ".npz": 
-  elif extension == ".h5":
+      return dictToNpz(fname,d)
-    return dictToH5(fname,d)
+    elif extension == ".h5":
-  else:
+      return dictToH5(fname,d)
-    raise ValueError("Extension must be h5 or npz") 
+    else:
+      raise ValueError("Extension must be h5 or npz") 
+  except Exception as e:
+    log.exception("Could not save %s"%fname)
 
 class DataStorage(dict):
   """ Storage for 1d integrated info """
-  def __init__(self,fileOrDict,default_name='pyfai_1d',default_ext='npz'):
+  def __init__(self,fileOrDict,recursive=True,
+               default_name='pyfai_1d',default_ext='npz'):
     if isinstance(fileOrDict,dict):
       self.filename = None
       d = fileOrDict
@@ -109,6 +113,11 @@ class DataStorage(dict):
       self.filename = fileOrDict
       d = read(fileOrDict)
 
+    if recursive:
+      for k in d.keys():
+        if not isinstance(d[k],DataStorage) and isinstance(d[k],dict):
+          d[k] = DataStorage(d[k])
+
     # allow accessing with .data, .delays, etc.
     for k,v in d.items(): setattr(self,k,v)
 
@@ -128,6 +137,31 @@ class DataStorage(dict):
     delattr(self,key)
     super().__delitem__(key)
 
+  def __str__(self):
+    keys = list(self.keys())
+    keys.sort()
+    return "DataStorage obj containing: %s" % ",".join(keys)
+
+  def __repr__(self):
+    keys = list(self.keys())
+    keys.sort()
+    nchars = max(map(len,keys))
+    fmt = "%%%ds %%s" % (nchars)
+    s = ["DataStorage obj containing (sorted): ",]
+    for k in keys:
+      if isinstance(self[k],np.ndarray):
+        value_str = "array %s"% "x".join(map(str,self[k].shape))
+      elif isinstance(self[k],DataStorage):
+        value_str = str(self[k])[:50] + "..."
+      elif isinstance(self[k],(str,DataStorage)):
+        value_str = self[k][:50] + "..."
+      elif self[k] is None:
+        value_str = "None"
+      else:
+        value_str = str(self[k])
+      s.append( fmt % (k,value_str) ) 
+    return "\n".join(s)
+
   def save(self,fname=None):
     if fname is None: fname = self.filename
     assert fname is not None

xray/utils.py

@@ -9,6 +9,7 @@ import os
 import glob
 import pathlib
 import re
+import numbers
 from . import storage as storage
 
 try:
@@ -72,6 +73,30 @@ def removeExt(fname):
 def getBasename(fname):
   return removeExt(os.path.basename(fname));
 
+def findSlice(array,lims):
+  start = np.ravel(np.argwhere(array>lims[0]))[0]
+  stop  = np.ravel(np.argwhere(array<lims[1]))[-1]
+  return slice(int(start),int(stop))
+
+def removeBackground(x,data,xlims=None,max_iter=100,background_regions=[],**kw):
+  from dualtree import dualtree
+  if data.ndim == 1: data = data[np.newaxis,:]
+  if xlims is not None:
+    idx = findSlice(x,xlims)
+    x = x[idx]
+    data = data[:,idx].copy()
+  else:
+    data = data.copy(); # create local copy
+  # has to be a list of lists ..
+  if background_regions != [] and isinstance(background_regions[0],numbers.Real):
+    background_regions = [background_regions,]
+  background_regions = [findSlice(x,brange) for brange in background_regions]
+  for i in range(len(data)):
+    data[i] = data[i] - dualtree.baseline(data[i],max_iter=max_iter,
+                        background_regions=background_regions,**kw)
+  return x,np.squeeze(data)
+  
+
 def plotdata(q,data,x=None,plot=True,showTrend=True,title=None,clim='auto'):
   if not (plot or showTrend): return
   if x is None: x = np.arange(data.shape[0])
@@ -108,7 +133,7 @@ def plotdiffs(q,diffs,t,select=None,err=None,absSignal=None,absSignalScale=10,
     indices = range(len(t))
   lines = []
   if absSignal is not None:
-    line = plt.plot(q,absSignal/absSignalScale,
+    line = plt.plot(q,absSignal/absSignalScale,lw=3,
                     color='k',label="absSignal/%s"%str(absSignalScale))[0]
     lines.append(line)
   for idiff in indices: