storage and xray modules live now on their own
This commit is contained in:
marco cammarata
parent
422c30e9b6
commit
5d14890e82
298
storage.py
298
storage.py
|
|
@ -1,298 +0,0 @@
|
||||||
""" npy/npz/hdf5 file based storage;
|
|
||||||
this modules adds the possibility to dump and load objects in files and
|
|
||||||
a more convenient was of accessing the data via the .attributedict thanks
|
|
||||||
to the DataStorage class """
|
|
||||||
import numpy as np
|
|
||||||
import os
|
|
||||||
import h5py
|
|
||||||
import collections
|
|
||||||
import logging
|
|
||||||
log = logging.getLogger(__name__)
|
|
||||||
|
|
||||||
_array_cache = {}
|
|
||||||
|
|
||||||
def unwrapArray(a,recursive=True,readH5pyDataset=True):
|
|
||||||
""" This function takes an object (like a dictionary) and recursively
|
|
||||||
unwraps it solving issues like:
|
|
||||||
* the fact that many objects are packaged as 0d array
|
|
||||||
This funciton has also some specific hack for handling h5py limits:
|
|
||||||
* handle the None python object
|
|
||||||
* numpy unicode ...
|
|
||||||
"""
|
|
||||||
# is h5py dataset convert to array
|
|
||||||
if isinstance(a,h5py.Dataset) and readH5pyDataset: a = a[...]
|
|
||||||
if isinstance(a,h5py.Dataset) and a.shape == (): a = a[...]
|
|
||||||
if isinstance(a,h5py.Group) and "IS_LIST_OF_ARRAYS" in a.attrs:
|
|
||||||
items = list(a.keys())
|
|
||||||
items.sort()
|
|
||||||
a = np.asarray( [a[item][...] for item in items] )
|
|
||||||
if isinstance(a,np.ndarray) and a.ndim == 0 : a = a.item()
|
|
||||||
if isinstance(a,np.ndarray) and a.dtype.char == "S": a = a.astype(str)
|
|
||||||
if recursive:
|
|
||||||
if "items" in dir(a): # dict, h5py groups, npz file
|
|
||||||
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,tuple)):
|
|
||||||
a = [unwrapArray(element) for element in a]
|
|
||||||
else:
|
|
||||||
pass
|
|
||||||
if isinstance(a,dict): a = DataStorage(a)
|
|
||||||
# restore None that cannot be saved in h5py
|
|
||||||
if isinstance(a,str) and a == "NONE_PYTHON_OBJECT": a = None
|
|
||||||
# h5py can't save numpy unicode
|
|
||||||
if isinstance(a,np.ndarray) and a.dtype.char == "S": a = a.astype(str)
|
|
||||||
return a
|
|
||||||
|
|
||||||
def dictToH5Group(d,group,link_copy=True):
|
|
||||||
""" helper function that transform (recursive) a dictionary into an
|
|
||||||
hdf group by creating subgroups
|
|
||||||
link_copy = True, tries to save space in the hdf file by creating an internal link.
|
|
||||||
the current implementation uses memory though ...
|
|
||||||
"""
|
|
||||||
global _array_cache
|
|
||||||
for key,value in d.items():
|
|
||||||
TOTRY = True
|
|
||||||
if isinstance(value,(list,tuple)): value = np.asarray(value)
|
|
||||||
if isinstance(value,dict):
|
|
||||||
group.create_group(key)
|
|
||||||
dictToH5Group(value,group[key],link_copy=link_copy)
|
|
||||||
elif value is None:
|
|
||||||
group[key] = "NONE_PYTHON_OBJECT"
|
|
||||||
elif isinstance(value,np.ndarray):
|
|
||||||
# take care of unicode (h5py can't handle numpy unicode arrays)
|
|
||||||
if value.dtype.char == "U": value = np.asarray([vv.encode('ascii') for vv in value])
|
|
||||||
# check if it is list of array
|
|
||||||
elif isinstance(value,np.ndarray) and value.ndim == 1 and isinstance(value[0],np.ndarray):
|
|
||||||
group.create_group(key)
|
|
||||||
group[key].attrs["IS_LIST_OF_ARRAYS"] = True
|
|
||||||
for index,array in enumerate(value): dictToH5Group( { "index%010d"%index : array},group[key],link_copy=link_copy );
|
|
||||||
TOTRY = False; # don't even try to save as generic call group[key]=value
|
|
||||||
if link_copy:
|
|
||||||
found_address = None
|
|
||||||
for address,(file_handle,array) in _array_cache.items():
|
|
||||||
if np.array_equal(array,value) and group.file == file_handle:
|
|
||||||
log.info("Found array in cache, asked for %s/%s, found as %s"%(group.name,key,address))
|
|
||||||
found_address = address
|
|
||||||
break
|
|
||||||
if found_address is not None:
|
|
||||||
value = group.file[found_address]
|
|
||||||
try:
|
|
||||||
if TOTRY:
|
|
||||||
group[key] = value
|
|
||||||
if link_copy:
|
|
||||||
log.info("Addind array %s to cache"%(group[key].name))
|
|
||||||
_array_cache[ group[key].name ] = (group.file,value)
|
|
||||||
except Exception as e:
|
|
||||||
log.warning("Can't save %s, error was %s"%(key,e))
|
|
||||||
# try saving everything else that is not dict or array
|
|
||||||
else:
|
|
||||||
try:
|
|
||||||
group[key] = value
|
|
||||||
except Exception as e:
|
|
||||||
log.error("Can't save %s, error was %s"%(key,e))
|
|
||||||
|
|
||||||
|
|
||||||
def dictToH5(h5,d,link_copy=False):
|
|
||||||
""" Save a dictionary into an hdf5 file
|
|
||||||
TODO: add capability of saving list of array
|
|
||||||
h5py is not capable of handling dictionaries natively"""
|
|
||||||
h5 = h5py.File(h5,mode="w")
|
|
||||||
# group = h5.create_group("/")
|
|
||||||
dictToH5Group(d,h5["/"],link_copy=link_copy)
|
|
||||||
h5.close()
|
|
||||||
|
|
||||||
def h5ToDict(h5,readH5pyDataset=True):
|
|
||||||
""" Read a hdf5 file into a dictionary """
|
|
||||||
with h5py.File(h5,"r") as h:
|
|
||||||
ret = unwrapArray(h,recursive=True,readH5pyDataset=readH5pyDataset)
|
|
||||||
return ret
|
|
||||||
|
|
||||||
def npzToDict(npzFile):
|
|
||||||
with np.load(npzFile) as npz: d = dict(npz)
|
|
||||||
d = unwrapArray(d,recursive=True)
|
|
||||||
return d
|
|
||||||
|
|
||||||
def npyToDict(npyFile):
|
|
||||||
d = unwrapArray( np.load(npyFile).item() ,recursive=True)
|
|
||||||
return d
|
|
||||||
|
|
||||||
def dictToNpz(npzFile,d): np.savez(npzFile,**d)
|
|
||||||
def dictToNpy(npyFile,d): np.save(npyFile,d)
|
|
||||||
|
|
||||||
def objToDict(o,recursive=True):
|
|
||||||
""" convert a DictWrap to a dictionary (useful for saving); it should work for other objects too
|
|
||||||
TODO: this function does not catch a list of DataStorage instances like
|
|
||||||
objToDict( ( DataStorage(), DataStorage() ) )
|
|
||||||
is not converted !!
|
|
||||||
"""
|
|
||||||
if "items" not in dir(o): return o
|
|
||||||
d = dict()
|
|
||||||
for k,v in o.items():
|
|
||||||
try:
|
|
||||||
d[k] = objToDict( v )
|
|
||||||
except Exception as e:
|
|
||||||
log.info("In objToDict, could not convert key %s to dict, error was"%\
|
|
||||||
(k,e))
|
|
||||||
d[k] = v
|
|
||||||
return d
|
|
||||||
|
|
||||||
|
|
||||||
def read(fname):
|
|
||||||
extension = os.path.splitext(fname)[1]
|
|
||||||
log.info("Reading storage file %s"%fname)
|
|
||||||
if extension == ".npz":
|
|
||||||
return DataStorage(npzToDict(fname))
|
|
||||||
elif extension == ".npy":
|
|
||||||
return DataStorage(npyToDict(fname))
|
|
||||||
elif extension == ".h5":
|
|
||||||
return DataStorage(h5ToDict(fname))
|
|
||||||
else:
|
|
||||||
raise ValueError("Extension must be h5, npy or npz, it was %s"%extension)
|
|
||||||
|
|
||||||
def save(fname,d,link_copy=True):
|
|
||||||
""" link_copy is used by hdf5 saving only, it allows to creat link of identical arrays (saving space) """
|
|
||||||
# make sure the object is dict (recursively) this allows reading it
|
|
||||||
# without the DataStorage module
|
|
||||||
d = objToDict(d,recursive=True)
|
|
||||||
d['filename'] = fname
|
|
||||||
extension = os.path.splitext(fname)[1]
|
|
||||||
log.info("Saving storage file %s"%fname)
|
|
||||||
try:
|
|
||||||
if extension == ".npz":
|
|
||||||
return dictToNpz(fname,d)
|
|
||||||
elif extension == ".h5":
|
|
||||||
return dictToH5(fname,d,link_copy=link_copy)
|
|
||||||
elif extension == ".npy":
|
|
||||||
return dictToNpy(fname,d)
|
|
||||||
else:
|
|
||||||
raise ValueError("Extension must be h5, npy or npz, it was %s"%extension)
|
|
||||||
except Exception as e:
|
|
||||||
log.exception("Could not save %s"%fname)
|
|
||||||
|
|
||||||
|
|
||||||
class DataStorage(dict):
|
|
||||||
""" Storage for dict like object.
|
|
||||||
recursive : bool
|
|
||||||
recursively convert dict-like objects to DataStorage
|
|
||||||
It can save data to file (format npy,npz or h5)
|
|
||||||
|
|
||||||
To initialize it:
|
|
||||||
|
|
||||||
data = DataStorage( a=(1,2,3),b="add",filename='store.npz' )
|
|
||||||
|
|
||||||
# recursively by default
|
|
||||||
# data.a will be a DataStorage instance
|
|
||||||
data = DataStorage( a=dict( b = 1)) );
|
|
||||||
|
|
||||||
# data.a will be a dictionary
|
|
||||||
data = DataStorage( a=dict( b = 1),recursive=False )
|
|
||||||
|
|
||||||
# reads from file if it exists
|
|
||||||
data = DataStorage( 'mysaveddata.npz' ) ;
|
|
||||||
|
|
||||||
DOES NOT READ FROM FILE (even if it exists)!!
|
|
||||||
data = DataStorage( filename = 'mysaveddata.npz' );
|
|
||||||
|
|
||||||
create empty storage (with default filename)
|
|
||||||
data = DataStorage()
|
|
||||||
"""
|
|
||||||
def __init__(self,*args,filename='data_storage.npz',recursive=True,**kwargs):
|
|
||||||
# self.filename = kwargs.pop('filename',"data_storage.npz")
|
|
||||||
self.filename = filename
|
|
||||||
self._recursive = recursive
|
|
||||||
# interpret kwargs as dict if there are
|
|
||||||
if len(kwargs) != 0:
|
|
||||||
fileOrDict = dict(kwargs)
|
|
||||||
elif len(kwargs)==0 and len(args)>0:
|
|
||||||
fileOrDict = args[0]
|
|
||||||
else:
|
|
||||||
fileOrDict = dict()
|
|
||||||
|
|
||||||
|
|
||||||
d = dict(); # data dictionary
|
|
||||||
if isinstance(fileOrDict,dict):
|
|
||||||
d = fileOrDict
|
|
||||||
elif isinstance(fileOrDict,str):
|
|
||||||
if os.path.isfile(fileOrDict):
|
|
||||||
d = read(fileOrDict)
|
|
||||||
else:
|
|
||||||
self.filename=fileOrDict
|
|
||||||
d = dict()
|
|
||||||
else:
|
|
||||||
raise ValueError("Invalid DataStorage definition")
|
|
||||||
|
|
||||||
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)
|
|
||||||
|
|
||||||
# allow accessing as proper dict
|
|
||||||
self.update( **dict(d) )
|
|
||||||
|
|
||||||
def __setitem__(self, key, value):
|
|
||||||
#print("__setitem__")
|
|
||||||
setattr(self,key,value)
|
|
||||||
super().__setitem__(key, value)
|
|
||||||
|
|
||||||
def __setattr__(self, key, value):
|
|
||||||
""" allows to add fields with data.test=4 """
|
|
||||||
# check if attr exists is essential (or it fails when defining an instance)
|
|
||||||
if hasattr(self,"_recursive") and self._recursive and \
|
|
||||||
isinstance(value,(dict,collections.OrderedDict)): value = DataStorage(value)
|
|
||||||
super().__setitem__(key, value)
|
|
||||||
super().__setattr__(key,value)
|
|
||||||
|
|
||||||
def __delitem__(self, key):
|
|
||||||
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()
|
|
||||||
if len(keys) == 0: return "Empty DataStorage"
|
|
||||||
nchars = max(map(len,keys))
|
|
||||||
fmt = "%%%ds %%s" % (nchars)
|
|
||||||
s = ["DataStorage obj containing (sorted): ",]
|
|
||||||
for k in keys:
|
|
||||||
if k[0] == "_": continue
|
|
||||||
obj = self[k]
|
|
||||||
if ( (isinstance(obj,np.ndarray) and obj.ndim == 1) or isinstance(obj,(list,tuple))) and all( [isinstance(v,np.ndarray) for v in obj]):
|
|
||||||
value_str = "list of arrays, shapes " + ",".join([str(v.shape) for v in obj[:5]]) + " ..."
|
|
||||||
elif isinstance(obj,np.ndarray):
|
|
||||||
value_str = "array, size %s, type %s"% ("x".join(map(str,obj.shape)),obj.dtype)
|
|
||||||
elif isinstance(obj,DataStorage):
|
|
||||||
value_str = str(obj)[:50]
|
|
||||||
elif isinstance(obj,(str,DataStorage)):
|
|
||||||
value_str = obj[:50]
|
|
||||||
elif self[k] is None:
|
|
||||||
value_str = "None"
|
|
||||||
else:
|
|
||||||
value_str = str(self[k])
|
|
||||||
if len(str(obj))>50: value_str += " ..."
|
|
||||||
s.append( fmt % (k,value_str) )
|
|
||||||
return "\n".join(s)
|
|
||||||
|
|
||||||
def keys(self):
|
|
||||||
keys = list(super().keys())
|
|
||||||
keys = [k for k in keys if k != 'filename' ]
|
|
||||||
keys = [k for k in keys if k[0] != '_' ]
|
|
||||||
return keys
|
|
||||||
|
|
||||||
def save(self,fname=None,link_copy=False):
|
|
||||||
""" link_copy: only works in hfd5 format
|
|
||||||
save space by creating link when identical arrays are found,
|
|
||||||
it slows down the saving (3 or 4 folds) but saves A LOT of space
|
|
||||||
when saving different dataset together (since it does not duplicate
|
|
||||||
internal pyfai matrices
|
|
||||||
"""
|
|
||||||
if fname is None: fname = self.filename
|
|
||||||
assert fname is not None
|
|
||||||
save(fname,self,link_copy=link_copy)
|
|
||||||
|
|
@ -1,24 +0,0 @@
|
||||||
# xray
|
|
||||||
|
|
||||||
Different utilities to work with time resolved solution/powder scattering data
|
|
||||||
They should be extended to absorption soon
|
|
||||||
|
|
||||||
The idea is an hierarchical structure of modules
|
|
||||||
|
|
||||||
> user → beamline → workhorses
|
|
||||||
|
|
||||||
## Workhorses
|
|
||||||
azav.py: uses pyFAI to convert images into 1d curves
|
|
||||||
dataReduction.py: works on the normalization, referencing and averaging.
|
|
||||||
|
|
||||||
Both modules uses some common utilities in utils.py including file based storage for persistency and convenience (numpy or hdf5 based);
|
|
||||||
Indeed the data_storage class wraps around dictionaries allowing accessing with .attribute syntax.
|
|
||||||
In utils plotfuncitons are presents
|
|
||||||
|
|
||||||
## beamline
|
|
||||||
this step is needed to read the diagnostic information (time delay, monitors, etc). An example for current id9 data collection macros is given (id9.py). Essentially it means creating a dictionary where each key will be added to the main data dictionary (that has keys like q, data, files, etc.)
|
|
||||||
|
|
||||||
## user
|
|
||||||
at the user level, a rather simple macro like the one provided as an example (example_main_tiox.py) should be sufficient for most needs.
|
|
||||||
|
|
||||||
|
|
||||||
|
|
@ -1,12 +0,0 @@
|
||||||
from .. import storage
|
|
||||||
from . import utils
|
|
||||||
from . import mask
|
|
||||||
from . import peaks
|
|
||||||
from . import cell
|
|
||||||
from . import filters
|
|
||||||
from . import id9
|
|
||||||
try:
|
|
||||||
from . import azav
|
|
||||||
except ImportError:
|
|
||||||
print("Can't import azav ...")
|
|
||||||
from . import dataReduction
|
|
||||||
465
xray/azav.py
465
xray/azav.py
|
|
@ -1,465 +0,0 @@
|
||||||
from __future__ import print_function,division
|
|
||||||
|
|
||||||
import logging
|
|
||||||
log = logging.getLogger(__name__)
|
|
||||||
|
|
||||||
|
|
||||||
import numpy as np
|
|
||||||
np.seterr(all='ignore')
|
|
||||||
import os
|
|
||||||
import collections
|
|
||||||
import glob
|
|
||||||
import pathlib
|
|
||||||
from . import storage
|
|
||||||
from . import utils
|
|
||||||
from . import filters
|
|
||||||
import re
|
|
||||||
import fabio
|
|
||||||
import pyFAI
|
|
||||||
|
|
||||||
try:
|
|
||||||
import matplotlib.pyplot as plt
|
|
||||||
except ImportError:
|
|
||||||
log.warn("Can't import matplotlib !")
|
|
||||||
|
|
||||||
def _read(fname):
|
|
||||||
""" read data from file using fabio """
|
|
||||||
f = fabio.open(fname)
|
|
||||||
data = f.data
|
|
||||||
del f; # close file
|
|
||||||
return data
|
|
||||||
|
|
||||||
def read(fnames):
|
|
||||||
""" read data from file(s) using fabio """
|
|
||||||
if isinstance(fnames,str):
|
|
||||||
data = _read(fnames)
|
|
||||||
else:
|
|
||||||
# read one image to know img size
|
|
||||||
temp = _read(fnames[0])
|
|
||||||
shape = [len(fnames),]+list(temp.shape)
|
|
||||||
data = np.empty(shape)
|
|
||||||
data[0] = temp
|
|
||||||
for i in range(1,len(fnames)): data[i] = _read(fnames[i])
|
|
||||||
return data
|
|
||||||
|
|
||||||
def ai_as_dict(ai):
|
|
||||||
""" ai is a pyFAI azimuthal intagrator"""
|
|
||||||
methods = dir(ai)
|
|
||||||
methods = [m for m in methods if m.find("get_") == 0]
|
|
||||||
names = [m[4:] for m in methods]
|
|
||||||
values = [getattr(ai,m)() for m in methods]
|
|
||||||
ret = dict( zip(names,values) )
|
|
||||||
ret["detector"] = ai.detector.get_name()
|
|
||||||
return ret
|
|
||||||
|
|
||||||
def ai_as_str(ai):
|
|
||||||
""" ai is a pyFAI azimuthal intagrator"""
|
|
||||||
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,dezinger=None):
|
|
||||||
""" ai is a pyFAI azimuthal intagrator
|
|
||||||
it can be defined with pyFAI.load(ponifile)
|
|
||||||
dezinger: None or float (used as percentile of ai.separate)
|
|
||||||
mask: True are points to be masked out """
|
|
||||||
# force float to be sure of type casting for img
|
|
||||||
if isinstance(dark,int): dark = float(dark);
|
|
||||||
if imgs.ndim == 2: imgs = (imgs,)
|
|
||||||
out_i = np.empty( ( len(imgs), npt_radial) )
|
|
||||||
out_s = np.empty( ( len(imgs), npt_radial) )
|
|
||||||
for _i,img in enumerate(imgs):
|
|
||||||
if dezinger is not None and dezinger > 0:
|
|
||||||
_,img=ai.separate(img,npt_rad=npt_radial,npt_azim=512,unit='q_A^-1',
|
|
||||||
method=method,mask=mask,percentile=dezinger)
|
|
||||||
q,i, sig = ai.integrate1d(img-dark, npt_radial, mask= mask, safe = safe,\
|
|
||||||
unit="q_A^-1", method = method, error_model = "poisson",
|
|
||||||
polarization_factor = polCorr)
|
|
||||||
out_i[_i] = i
|
|
||||||
out_s[_i] = sig
|
|
||||||
return q,np.squeeze(out_i),np.squeeze(out_s)
|
|
||||||
|
|
||||||
def do2d(ai, imgs, mask = None, npt_radial = 600, npt_azim=360,method = 'csr',safe=True,dark=10., polCorr = 1):
|
|
||||||
""" ai is a pyFAI azimuthal intagrator
|
|
||||||
it can be defined with pyFAI.load(ponifile)
|
|
||||||
mask: True are points to be masked out """
|
|
||||||
# force float to be sure of type casting for img
|
|
||||||
if isinstance(dark,int): dark = float(dark);
|
|
||||||
if imgs.ndim == 2: imgs = (imgs,)
|
|
||||||
out = np.empty( ( len(imgs), npt_azim,npt_radial) )
|
|
||||||
for _i,img in enumerate(imgs):
|
|
||||||
i2d,q,azTheta = ai.integrate2d(img-dark, npt_radial, npt_azim=npt_azim,
|
|
||||||
mask= mask, safe = safe,unit="q_A^-1", method = method,
|
|
||||||
polarization_factor = polCorr )
|
|
||||||
out[_i] = i2d
|
|
||||||
return q,azTheta,np.squeeze(out)
|
|
||||||
|
|
||||||
def getAI(poni=None,folder=None,**kwargs):
|
|
||||||
""" get AzimuthalIntegrator instance:
|
|
||||||
→ if poni is a string, it is used as filename to read.
|
|
||||||
in this case if folder is given it is used (together with all its
|
|
||||||
subfolder) as search path (along with ./ and home folder)
|
|
||||||
→ kwargs if present can be used to define (or override) parameters from files,
|
|
||||||
dist,xcen,ycen,poni1,poni2,rot1,rot2,rot3,pixel1,pixel2,splineFile,
|
|
||||||
detector,wavelength
|
|
||||||
"""
|
|
||||||
if isinstance(poni,pyFAI.azimuthalIntegrator.AzimuthalIntegrator):
|
|
||||||
ai = poni
|
|
||||||
elif isinstance(poni,str):
|
|
||||||
# look is file exists in cwd
|
|
||||||
if os.path.isfile(poni):
|
|
||||||
fname = poni
|
|
||||||
# if file does not exist look for one with that name around
|
|
||||||
else:
|
|
||||||
# build search paths
|
|
||||||
folders = []
|
|
||||||
if folder is not None:
|
|
||||||
temp = os.path.abspath(folder)
|
|
||||||
path = pathlib.Path(temp)
|
|
||||||
folders = [ str(path), ]
|
|
||||||
for p in path.parents: folders.append(str(p))
|
|
||||||
folders.append( "./" )
|
|
||||||
folders.append( os.path.expanduser("~/") )
|
|
||||||
# look for file
|
|
||||||
for path in folders:
|
|
||||||
fname = path + "/" + poni
|
|
||||||
if os.path.isfile(fname):
|
|
||||||
log.info("Found poni file %s",fname)
|
|
||||||
break
|
|
||||||
else:
|
|
||||||
log.debug("Could not poni file %s",fname)
|
|
||||||
ai = pyFAI.load(fname)
|
|
||||||
else:
|
|
||||||
ai = pyFAI.azimuthalIntegrator.AzimuthalIntegrator()
|
|
||||||
for par,value in kwargs.items(): setattr(ai,par,value)
|
|
||||||
# provide xcen and ycen for convenience (note: xcen changes poni2
|
|
||||||
# and ycen changes poni1)
|
|
||||||
if 'xcen' in kwargs: ai.poni2 = kwargs['xcen'] * ai.pixel2
|
|
||||||
if 'ycen' in kwargs: ai.poni1 = kwargs['ycen'] * ai.pixel1
|
|
||||||
ai.reset(); # needed in case of overridden parameters
|
|
||||||
return ai
|
|
||||||
|
|
||||||
g_mask_str = re.compile("(\w)\s*(<|>)\s*(\d+)")
|
|
||||||
|
|
||||||
def _interpretMask(mask,shape=None):
|
|
||||||
"""
|
|
||||||
if mask is an existing filename, returns it
|
|
||||||
if mask is a string like [x|y] [<|>] int;
|
|
||||||
for example y>500 will dis-regard out for y>500
|
|
||||||
"""
|
|
||||||
maskout = None
|
|
||||||
if isinstance(mask,str) and os.path.isfile(mask):
|
|
||||||
maskout = read(mask).astype(np.bool)
|
|
||||||
elif isinstance(mask,str) and not os.path.isfile(mask):
|
|
||||||
err_msg = ValueError("The string '%s' could not be interpreted as simple\
|
|
||||||
mask; it should be something like x>10"%mask)
|
|
||||||
assert shape is not None
|
|
||||||
# interpret string
|
|
||||||
maskout = np.zeros(shape,dtype=bool)
|
|
||||||
match = g_mask_str.match(mask)
|
|
||||||
if match is None: raise err_msg
|
|
||||||
(axis,sign,lim) = match.groups()
|
|
||||||
if axis not in ("x","y"): raise err_msg
|
|
||||||
if sign not in (">","<"): raise err_msg
|
|
||||||
lim = int(lim)
|
|
||||||
idx = slice(lim,None) if sign == ">" else slice(None,lim)
|
|
||||||
if axis == 'y':
|
|
||||||
maskout[idx,:] = True
|
|
||||||
else:
|
|
||||||
maskout[:,idx] = True
|
|
||||||
elif isinstance(mask,np.ndarray):
|
|
||||||
maskout = mask.astype(np.bool)
|
|
||||||
elif mask is None:
|
|
||||||
assert shape is not None
|
|
||||||
maskout = np.zeros(shape,dtype=bool)
|
|
||||||
else:
|
|
||||||
maskout = None
|
|
||||||
raise ValueError("Could not interpret %s as mask input"%mask)
|
|
||||||
|
|
||||||
if shape is not None and maskout.shape != shape:
|
|
||||||
raise ValueError("The mask shape %s does not match the shape given as\
|
|
||||||
argument %s"%(maskout.shape,shape))
|
|
||||||
return maskout
|
|
||||||
|
|
||||||
|
|
||||||
def interpretMask(masks,shape=None):
|
|
||||||
"""
|
|
||||||
if masks is a list of masks, eachone can be:
|
|
||||||
* an existing filename
|
|
||||||
* a string like [x|y] [<|>] int;
|
|
||||||
"""
|
|
||||||
if not isinstance( masks, (list,tuple,np.ndarray) ):
|
|
||||||
masks = (masks,)
|
|
||||||
masks = [_interpretMask(mask,shape) for mask in masks]
|
|
||||||
# put them all together
|
|
||||||
mask = masks[0]
|
|
||||||
for m in masks[1:]:
|
|
||||||
mask = np.logical_or(mask,m)
|
|
||||||
return mask
|
|
||||||
|
|
||||||
def removeBackground(data,qlims=(0,10),max_iter=30,background_regions=[],force=False,
|
|
||||||
storageFile=None,save=True,**removeBkg):
|
|
||||||
""" similar function to the zray.utils one, this works on dataset created by
|
|
||||||
doFolder """
|
|
||||||
idx = utils.findSlice(data.orig.q,qlims)
|
|
||||||
# see if there are some to do ...
|
|
||||||
if force:
|
|
||||||
idx_start = 0
|
|
||||||
else:
|
|
||||||
idx_start = len(data.data)
|
|
||||||
if idx_start < len(data.orig.data):
|
|
||||||
_q,_data = utils.removeBackground(data.orig.q[idx],data.orig.data[:,idx],
|
|
||||||
max_iter=max_iter,background_regions=background_regions,**removeBkg)
|
|
||||||
data.q = _q
|
|
||||||
data.data = np.concatenate( (data.data,_data ) )
|
|
||||||
data.err = np.concatenate( (data.err ,data.err[idx_start,idx] ) )
|
|
||||||
if save: data.save(storageFile); # if None uses .filename
|
|
||||||
return data
|
|
||||||
|
|
||||||
|
|
||||||
def doFolder(folder,files='*.edf*',nQ = 1500,force=False,mask=None,dark=10,
|
|
||||||
norm='auto',save_pyfai=False,saveChi=True,poni='pyfai.poni',
|
|
||||||
storageFile='auto',save=True,logDict=None,dezinger=None,skip_first=0,last=None):
|
|
||||||
""" calc 1D curves from files in folder, returning a dictionary of stuff
|
|
||||||
nQ : number of Q-points (equispaced)
|
|
||||||
force : if True, redo from beginning even if previous data are found
|
|
||||||
if False, do only new files
|
|
||||||
mask : can be a list of [filenames|array of booleans|mask string]
|
|
||||||
pixels that are True are dis-regarded
|
|
||||||
saveChi: self-explanatory
|
|
||||||
dezinger: None or 0 to disable; good value is ~50. Needs good center and mask
|
|
||||||
logDict: dictionary(-like) structure. has to have 'file' key
|
|
||||||
save_pyfai: store all pyfai's internal arrays (~110 MB)
|
|
||||||
poni : could be:
|
|
||||||
→ an AzimuthalIntegrator instance
|
|
||||||
→ a filename that will be look for in
|
|
||||||
1 'folder' first
|
|
||||||
2 in ../folder
|
|
||||||
3 in ../../folder
|
|
||||||
....
|
|
||||||
n-1 in pwd
|
|
||||||
n in homefolder
|
|
||||||
→ a dictionary (use to bootstrap an AzimuthalIntegrator using
|
|
||||||
AzimuthalIntegrator(**poni)
|
|
||||||
"""
|
|
||||||
if storageFile == 'auto': storageFile = folder + "/" + "pyfai_1d.h5"
|
|
||||||
|
|
||||||
if os.path.isfile(storageFile) and not force:
|
|
||||||
saved = storage.DataStorage(storageFile)
|
|
||||||
log.info("Found %d images in storage file"%saved.data.shape[0])
|
|
||||||
else:
|
|
||||||
saved = None
|
|
||||||
|
|
||||||
files = utils.getFiles(folder,files)
|
|
||||||
if logDict is not None:
|
|
||||||
files = [f for f in files if utils.getBasename(f) in logDict['file'] ]
|
|
||||||
files = files[skip_first:last]
|
|
||||||
|
|
||||||
|
|
||||||
if saved is not None:
|
|
||||||
files = [f for f in files if f not in saved["files"]]
|
|
||||||
log.info("Will do azimuthal integration for %d files"%(len(files)))
|
|
||||||
|
|
||||||
files = np.asarray(files)
|
|
||||||
basenames = np.asarray( [ utils.getBasename(file) for file in files] )
|
|
||||||
|
|
||||||
if len(files) > 0:
|
|
||||||
# which poni file to use:
|
|
||||||
ai = getAI(poni,folder)
|
|
||||||
|
|
||||||
|
|
||||||
shape = read(files[0]).shape
|
|
||||||
mask = interpretMask(mask,shape)
|
|
||||||
|
|
||||||
data = np.empty( (len(files),nQ) )
|
|
||||||
err = np.empty( (len(files),nQ) )
|
|
||||||
for ifname,fname in enumerate(files):
|
|
||||||
img = read(fname)
|
|
||||||
q,i,e = do1d(ai,img,mask=mask,npt_radial=nQ,dark=dark,dezinger=dezinger)
|
|
||||||
data[ifname] = i
|
|
||||||
err[ifname] = e
|
|
||||||
if saveChi:
|
|
||||||
chi_fname = utils.removeExt(fname) + ".chi"
|
|
||||||
utils.saveTxt(chi_fname,q,np.vstack((i,e)),info=ai_as_str(ai),overwrite=True)
|
|
||||||
|
|
||||||
if saved is not None:
|
|
||||||
files = np.concatenate( (saved["files"] ,basenames ) )
|
|
||||||
data = np.concatenate( (saved["data"] ,data ) )
|
|
||||||
err = np.concatenate( (saved["err"] ,err ) )
|
|
||||||
theta_rad = utils.qToTheta(q,wavelength=ai.wavelength)
|
|
||||||
theta_deg = utils.qToTheta(q,wavelength=ai.wavelength,asDeg=True)
|
|
||||||
orig = dict(data=data.copy(),err=err.copy(),q=q.copy())
|
|
||||||
ret = dict(q=q,folder=folder,files=files,data=data,err=err,
|
|
||||||
orig = orig, theta_rad = theta_rad, theta_deg=theta_deg,
|
|
||||||
pyfai=ai_as_dict(ai),pyfai_info=ai_as_str(ai),mask=mask)
|
|
||||||
if not save_pyfai:
|
|
||||||
ret['pyfai']['chia'] = None
|
|
||||||
ret['pyfai']['dssa'] = None
|
|
||||||
ret['pyfai']['q'] = None
|
|
||||||
ret['pyfai']['ttha'] = None
|
|
||||||
|
|
||||||
|
|
||||||
ret = storage.DataStorage(ret)
|
|
||||||
|
|
||||||
# add info from logDict if provided
|
|
||||||
if logDict is not None:
|
|
||||||
ret['log']=logDict
|
|
||||||
# sometime saving is not necessary (if one has to do it after subtracting background
|
|
||||||
if storageFile is not None and save: ret.save(storageFile)
|
|
||||||
else:
|
|
||||||
ret = saved
|
|
||||||
return ret
|
|
||||||
|
|
||||||
|
|
||||||
def removeBackground(data,qlims=(0,10),max_iter=30,background_regions=[],
|
|
||||||
storageFile=None,save=True,**removeBkg):
|
|
||||||
""" similar function to the zray.utils one, this works on dataset created by
|
|
||||||
doFolder """
|
|
||||||
idx = utils.findSlice(data.q,qlims)
|
|
||||||
data.q,data.data = utils.removeBackground(data.q[idx],data.data[:,idx],
|
|
||||||
max_iter=max_iter,background_regions=background_regions,**removeBkg)
|
|
||||||
data.err = data.err[idx]
|
|
||||||
if save: data.save(storageFile); # if None uses .filename
|
|
||||||
return data
|
|
||||||
|
|
||||||
def _calc_R(x,y, xc, yc):
|
|
||||||
""" calculate the distance of each 2D points from the center (xc, yc) """
|
|
||||||
return np.sqrt((x-xc)**2 + (y-yc)**2)
|
|
||||||
|
|
||||||
def _chi2(c, x, y):
|
|
||||||
""" calculate the algebraic distance between the data points and the mean
|
|
||||||
circle centered at c=(xc, yc) """
|
|
||||||
Ri = _calc_R(x, y, *c)
|
|
||||||
return Ri - Ri.mean()
|
|
||||||
|
|
||||||
def leastsq_circle(x,y):
|
|
||||||
from scipy import optimize
|
|
||||||
# coordinates of the barycenter
|
|
||||||
center_estimate = np.nanmean(x), np.nanmean(y)
|
|
||||||
center, ier = optimize.leastsq(_chi2, center_estimate, args=(x,y))
|
|
||||||
xc, yc = center
|
|
||||||
Ri = _calc_R(x, y, *center)
|
|
||||||
R = Ri.mean()
|
|
||||||
residu = np.sum((Ri - R)**2)
|
|
||||||
return xc, yc, R
|
|
||||||
|
|
||||||
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)
|
|
||||||
ai = pyFAI.azimuthalIntegrator.AzimuthalIntegrator(**kw)
|
|
||||||
fig_img,ax_img = plt.subplots(1,1)
|
|
||||||
fig_pyfai,ax_pyfai = plt.subplots(1,1)
|
|
||||||
fig_pyfai = plt.figure(2)
|
|
||||||
temp= ax_img.imshow(img)
|
|
||||||
plt.sca(ax_img); # set figure to use for mouse interaction
|
|
||||||
temp.set_clim( *np.percentile(img,(2,95) ) )
|
|
||||||
ans = ""
|
|
||||||
print("Enter 'end' when done")
|
|
||||||
while ans != "end":
|
|
||||||
if center is None:
|
|
||||||
print("Click on beam center:")
|
|
||||||
plt.sca(ax_img); # set figure to use for mouse interaction
|
|
||||||
center = plt.ginput()[0]
|
|
||||||
print("Selected center:",center)
|
|
||||||
ai.set_poni1(center[1]*psize)
|
|
||||||
ai.set_poni2(center[0]*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(center))
|
|
||||||
if reference is not None: ax_pyfai.axvline(reference)
|
|
||||||
plt.pause(0.01)
|
|
||||||
plt.draw()
|
|
||||||
plt.draw_all()
|
|
||||||
ans=input("Enter to continue with clinking or enter xc,yc values ")
|
|
||||||
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,
|
|
||||||
showTrend=False):
|
|
||||||
data = storage.DataStorage(fileOrFolder)
|
|
||||||
if isinstance(delays,slice):
|
|
||||||
idx = np.arange(data.delays.shape[0])[delays]
|
|
||||||
elif isinstance(delays,(int,float)):
|
|
||||||
idx = data.delays == float(delays)
|
|
||||||
else:
|
|
||||||
idx = data.delays < 0
|
|
||||||
if idx.sum() == 0:
|
|
||||||
print("No data with the current filter")
|
|
||||||
return None
|
|
||||||
i = data.data[idx]
|
|
||||||
q = data.q
|
|
||||||
if isinstance(norm,(tuple,list)):
|
|
||||||
idx = ( q>norm[0] ) & (q<norm[1])
|
|
||||||
norm = np.nanmean(i[:,idx],axis=1)
|
|
||||||
i = i/norm[:,np.newaxis]
|
|
||||||
if isinstance(norm,np.ndarray):
|
|
||||||
i = i/norm[:,np.newaxis]
|
|
||||||
title = "%s %s" % (fileOrFolder,str(delays))
|
|
||||||
utils.plotdata(q,i*scale,showTrend=showTrend,plot=plot,title=title)
|
|
||||||
if returnAll:
|
|
||||||
return q,i.mean(axis=0)*scale,i
|
|
||||||
else:
|
|
||||||
return q,i.mean(axis=0)*scale
|
|
||||||
|
|
||||||
#### Utilities for chi files ####
|
|
||||||
def chiRead(fname,scale=1):
|
|
||||||
q,i = np.loadtxt(fname,unpack=True,usecols=(0,1))
|
|
||||||
return q,i*scale
|
|
||||||
|
|
||||||
def chiPlot(fname,useTheta=False,E=12.4):
|
|
||||||
q,i = chiRead(fname)
|
|
||||||
lam = 12.4/E
|
|
||||||
theta = 2*180/3.14*np.arcsin(q*lam/4/3.14)
|
|
||||||
x = theta if useTheta else q
|
|
||||||
plt.plot(x,i,label=fname)
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
def chiAverage(folder,basename="",scale=1,norm=None,returnAll=False,plot=False,showTrend=False,clim='auto'):
|
|
||||||
files = glob.glob("%s/%s*chi"%(folder,basename))
|
|
||||||
files.sort()
|
|
||||||
print(files)
|
|
||||||
if len(files) == 0:
|
|
||||||
print("No file found (basename %s)" % basename)
|
|
||||||
return None
|
|
||||||
q,_ = chiRead(files[0])
|
|
||||||
i = np.asarray( [ chiRead(f)[1] for f in files ] )
|
|
||||||
if isinstance(norm,(tuple,list)):
|
|
||||||
idx = ( q>norm[0] ) & (q<norm[1])
|
|
||||||
norm = np.nanmean(i[:,idx],axis=1)
|
|
||||||
i = i/norm[:,np.newaxis]
|
|
||||||
elif isinstance(norm,np.ndarray):
|
|
||||||
i = i/norm[:,np.newaxis]
|
|
||||||
title = "%s %s" % (folder,basename)
|
|
||||||
utils.plotdata(q,i,plot=plot,showTrend=showTrend,title=title,clim=clim)
|
|
||||||
if (showTrend and plot): plt.subplot(1,2,1)
|
|
||||||
if showTrend:
|
|
||||||
plt.pcolormesh(np.arange(i.shape[0]),q,i.T)
|
|
||||||
plt.xlabel("image number, 0 being older")
|
|
||||||
plt.ylabel(r"q ($\AA^{-1}$)")
|
|
||||||
if (showTrend and plot): plt.subplot(1,2,2)
|
|
||||||
if plot:
|
|
||||||
plt.plot(q,i.mean(axis=0)*scale)
|
|
||||||
if (plot or showTrend):
|
|
||||||
plt.title(folder+"/"+basename)
|
|
||||||
if returnAll:
|
|
||||||
return q,i.mean(axis=0)*scale,i
|
|
||||||
else:
|
|
||||||
return q,i.mean(axis=0)*scale
|
|
||||||
98
xray/cell.py
98
xray/cell.py
|
|
@ -1,98 +0,0 @@
|
||||||
from __future__ import division,print_function
|
|
||||||
import collections
|
|
||||||
import itertools
|
|
||||||
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)
|
|
||||||
|
|
||||||
def reflection_list(self,maxQ=3,lim=10):
|
|
||||||
ret=dict()
|
|
||||||
# prepare hkl
|
|
||||||
i = range(-lim,lim+1)
|
|
||||||
prod = itertools.product( i,i,i )
|
|
||||||
hkl = np.asarray( list( itertools.product( i,i,i ) ) )
|
|
||||||
h,k,l = hkl.T
|
|
||||||
q = self.Q(h,k,l)
|
|
||||||
|
|
||||||
idx = q<maxQ;
|
|
||||||
q = q[idx]
|
|
||||||
hkl = hkl[idx]
|
|
||||||
|
|
||||||
qunique = np.unique(q)
|
|
||||||
ret = []
|
|
||||||
for qi in qunique:
|
|
||||||
reflec = hkl[ q == qi ]
|
|
||||||
ret.append( (qi,tuple(np.abs(reflec)[0]),len(reflec),reflec) )
|
|
||||||
return qunique,ret
|
|
||||||
|
|
||||||
# for h in range(-lim,lim+1):
|
|
||||||
# for j in range(-lim,lim+1):
|
|
||||||
|
|
||||||
|
|
||||||
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.):
|
|
||||||
Triclinic.__init__(self,a=a,b=b,c=c,alpha=90,beta=beta,gamma=90)
|
|
||||||
|
|
||||||
|
|
||||||
def plotReflections(cell_instance,ax=None,line_kw=dict(),text_kw=dict()):
|
|
||||||
import matplotlib.pyplot as plt
|
|
||||||
from matplotlib import lines
|
|
||||||
import matplotlib.transforms as transforms
|
|
||||||
_,refl_info = cell_instance.reflection_list()
|
|
||||||
if ax is None: ax = plt.gca()
|
|
||||||
|
|
||||||
# the x coords of this transformation are data, and the
|
|
||||||
# y coord are axes
|
|
||||||
trans = transforms.blended_transform_factory(ax.transData, ax.transAxes)
|
|
||||||
txt_kw = dict( horizontalalignment='center', rotation=45)
|
|
||||||
txt_kw.update(**text_kw)
|
|
||||||
for reflection in refl_info[1:]:
|
|
||||||
q,hkl,n,_ = reflection
|
|
||||||
line = lines.Line2D( [q,q],[1,1.1],transform=trans,**line_kw)
|
|
||||||
line.set_clip_on(False)
|
|
||||||
ax.add_line(line)
|
|
||||||
ax.text(q,1.15,str(hkl),transform=trans,**txt_kw)
|
|
||||||
|
|
||||||
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)
|
|
||||||
|
|
@ -1,261 +0,0 @@
|
||||||
from __future__ import print_function,division
|
|
||||||
|
|
||||||
import logging
|
|
||||||
log = logging.getLogger(__name__)
|
|
||||||
|
|
||||||
import numpy as np
|
|
||||||
np.seterr(all='ignore')
|
|
||||||
from . import utils
|
|
||||||
from . import storage
|
|
||||||
import os
|
|
||||||
|
|
||||||
def subtractReferences(i,idx_ref, useRatio = False):
|
|
||||||
""" given data in i (first index is shot num) and the indeces of the
|
|
||||||
references (idx_ref, array of integers) it interpolates the closest
|
|
||||||
reference data for each shot and subtracts it (or divides it, depending
|
|
||||||
on useRatio = [True|False];
|
|
||||||
Note: it works in place (i.e. it modifies i) """
|
|
||||||
iref=np.empty_like(i)
|
|
||||||
idx_ref = np.squeeze(idx_ref)
|
|
||||||
idx_ref = np.atleast_1d(idx_ref)
|
|
||||||
# sometime there is just one reference (e.g. sample scans)
|
|
||||||
if idx_ref.shape[0] == 1:
|
|
||||||
if useRatio:
|
|
||||||
return i/i[idx_ref]
|
|
||||||
else:
|
|
||||||
return i-i[idx_ref]
|
|
||||||
# references before first ref are "first ref"
|
|
||||||
iref[:idx_ref[0]] = i[idx_ref[0]]
|
|
||||||
# references after last ref are "last ref"
|
|
||||||
iref[idx_ref[-1]:] = i[idx_ref[-1]]
|
|
||||||
_ref = 0
|
|
||||||
for _i in range(idx_ref[0],idx_ref[-1]):
|
|
||||||
if _i in idx_ref: continue
|
|
||||||
idx_ref_before = idx_ref[_ref]
|
|
||||||
idx_ref_after = idx_ref[_ref+1]
|
|
||||||
ref_before = i[idx_ref_before]
|
|
||||||
ref_after = i[idx_ref_after]
|
|
||||||
weight_before = float(_i-idx_ref_before)/(idx_ref_after-idx_ref_before)
|
|
||||||
weight_after = 1-weight_before
|
|
||||||
# normal reference for an on chi, the weighted average
|
|
||||||
iref[_i] = weight_before*ref_before + weight_after*ref_after
|
|
||||||
if _i>=idx_ref_after-1: _ref += 1
|
|
||||||
log.debug("For image %d : %d-%d"%(_i,idx_ref_before,idx_ref_after))
|
|
||||||
# take care of the reference for the references ...
|
|
||||||
if len(idx_ref) > 2:
|
|
||||||
iref[idx_ref[0]] = i[idx_ref[1]]
|
|
||||||
iref[idx_ref[-1]] = i[idx_ref[-2]]
|
|
||||||
for _i in range(1,len(idx_ref)-1):
|
|
||||||
idx_ref_before = idx_ref[_i-1]
|
|
||||||
idx_ref_after = idx_ref[_i+1]
|
|
||||||
ref_before = i[idx_ref_before]
|
|
||||||
ref_after = i[idx_ref_after]
|
|
||||||
weight_before = float(idx_ref[_i]-idx_ref_before)/(idx_ref_after-idx_ref_before)
|
|
||||||
weight_after = 1-weight_before
|
|
||||||
# normal reference for an on chi, the weighted average
|
|
||||||
iref[idx_ref[_i]] = weight_before*ref_before + weight_after*ref_after
|
|
||||||
log.debug("For reference image %d : %d-%d"%(idx_ref[_i],idx_ref_before,idx_ref_after))
|
|
||||||
else:
|
|
||||||
#print(idx_ref)
|
|
||||||
#print(iref[idx_ref])
|
|
||||||
iref[idx_ref]=i[idx_ref[0]]
|
|
||||||
#print(iref[idx_ref])
|
|
||||||
if useRatio:
|
|
||||||
i /= iref
|
|
||||||
else:
|
|
||||||
i -= iref
|
|
||||||
return i
|
|
||||||
|
|
||||||
def averageScanPoints(scan,data,errAbs=None,isRef=None,lpower=None,useRatio=False,\
|
|
||||||
funcForAveraging=np.nanmean):
|
|
||||||
""" given scanpoints in 'scan' and corresponding data in 'data'
|
|
||||||
average all data corresponding the exactly the same scanpoint.
|
|
||||||
If the values in scan are coming from a readback, rounding might be
|
|
||||||
necessary.
|
|
||||||
No normalization is done inside this function
|
|
||||||
if isRef is provided must be a boolean array of the same shape as 'scan'
|
|
||||||
is there is at least one scanpoint marked as True, the data are
|
|
||||||
subtracted/divided by the interpolated reference
|
|
||||||
if lpower is provided the data is divided by it (how it is done depends
|
|
||||||
if one uses the ratio or not
|
|
||||||
funcForAveraging: is usually np.nanmean or np.nanmedian. it can be any
|
|
||||||
function that support axis=0 as keyword argument
|
|
||||||
"""
|
|
||||||
data = data.astype(np.float)
|
|
||||||
avData = np.nanmedian( data , axis = 0 )
|
|
||||||
|
|
||||||
if isRef is None: isRef = np.zeros( data.shape[0], dtype=bool )
|
|
||||||
assert data.shape[0] == isRef.shape[0]
|
|
||||||
|
|
||||||
# subtract reference only is there is at least one
|
|
||||||
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:
|
|
||||||
lpower = utils.reshapeToBroadcast(lpower,data)
|
|
||||||
if useRatio is False:
|
|
||||||
diff /= lpower
|
|
||||||
else:
|
|
||||||
diff = (data-1)/lpower+1
|
|
||||||
|
|
||||||
scan_pos = np.unique(scan)
|
|
||||||
shape_out = [len(scan_pos),] + list(diff.shape[1:])
|
|
||||||
ret = np.empty(shape_out)
|
|
||||||
err = np.empty(shape_out)
|
|
||||||
data_abs = np.empty(shape_out)
|
|
||||||
diffsInScanPoint = []
|
|
||||||
chi2_0 = []
|
|
||||||
for i,t in enumerate(scan_pos):
|
|
||||||
shot_idx = (scan == t)
|
|
||||||
|
|
||||||
# select data for the scan point
|
|
||||||
diff_for_scan = diff[shot_idx]
|
|
||||||
#if errAbs is not None:
|
|
||||||
# noise = np.nanmean(errAbs[shot_idx],axis = 0)
|
|
||||||
#else:
|
|
||||||
noise = np.nanstd(diff_for_scan, axis = 0)
|
|
||||||
|
|
||||||
# if it is the reference take only every second ...
|
|
||||||
if np.all( shot_idx == isRef ):
|
|
||||||
diff_for_scan = diff_for_scan[::2]
|
|
||||||
|
|
||||||
diffsInScanPoint.append( diff_for_scan )
|
|
||||||
|
|
||||||
# calculate average
|
|
||||||
ret[i] = funcForAveraging(diff_for_scan,axis=0)
|
|
||||||
data_abs[i] = funcForAveraging(data[shot_idx],axis=0)
|
|
||||||
|
|
||||||
# calculate chi2 of different repetitions
|
|
||||||
chi2 = np.power( (diff_for_scan - ret[i])/noise,2)
|
|
||||||
# sum over all axis but first
|
|
||||||
for _ in range(diff_for_scan.ndim-1):
|
|
||||||
chi2 = np.nansum( chi2, axis=-1 )
|
|
||||||
|
|
||||||
# store chi2_0
|
|
||||||
chi2_0.append( chi2/ret[i].size )
|
|
||||||
|
|
||||||
# store error of mean
|
|
||||||
err[i] = noise/np.sqrt(shot_idx.sum())
|
|
||||||
ret = dict(scan=scan_pos,data=ret,dataUnmasked=ret.copy(),err=err,
|
|
||||||
errUnmasked=err.copy(),chi2_0=chi2_0,diffsInScanPoint=diffsInScanPoint,
|
|
||||||
dataAbsAvNeg = avNeg, dataAsAbs=ret+avNeg,errAbs=errAbs,
|
|
||||||
dataAbsAvAll=avData,dataAbsAvScanPoint=data_abs,dataAbs=data.copy())
|
|
||||||
ret = storage.DataStorage(ret)
|
|
||||||
return ret
|
|
||||||
|
|
||||||
|
|
||||||
def calcTimeResolvedSignal(scan,data,err=None,reference="min",monitor=None,q=None,
|
|
||||||
saveTxt=True,folder="./",**kw):
|
|
||||||
"""
|
|
||||||
reference: can be 'min', 'max', a float|integer or an array of booleans
|
|
||||||
q : is needed if monitor is a tuple|list
|
|
||||||
monitor : normalization vector (if it is interpreted a list it is
|
|
||||||
interpreted as q-range normalization)
|
|
||||||
saveTxt : will save txt outputfiles (diff_av_*) in folder
|
|
||||||
other keywords are passed to averageScanPoints
|
|
||||||
"""
|
|
||||||
if reference == "min":
|
|
||||||
isRef = (scan == scan.min())
|
|
||||||
elif reference == "max":
|
|
||||||
isRef = (scan == scan.max())
|
|
||||||
elif isinstance(reference,(float,int)):
|
|
||||||
isRef = (scan == reference)
|
|
||||||
else:
|
|
||||||
isRef = reference
|
|
||||||
# normalize if needed
|
|
||||||
if monitor is not None:
|
|
||||||
if isinstance(monitor,(tuple,list)):
|
|
||||||
assert q is not None
|
|
||||||
assert data.ndim == 2
|
|
||||||
idx = (q>= monitor[0]) & (q<= monitor[1])
|
|
||||||
monitor = np.nanmedian(data[:,idx],axis=1)
|
|
||||||
monitor = utils.reshapeToBroadcast(monitor,data)
|
|
||||||
data = data/monitor
|
|
||||||
if err is not None: err = err/monitor
|
|
||||||
ret = averageScanPoints(scan,data,errAbs=err,isRef=isRef,**kw)
|
|
||||||
if q is not None: ret["q"] = q
|
|
||||||
return ret
|
|
||||||
|
|
||||||
def saveTxt(folder,data,delayToStr=True,basename='auto',info="",**kw):
|
|
||||||
""" data must be a DataStorage instance """
|
|
||||||
# folder ends usually with sample/run so use the last two subfolders
|
|
||||||
# the abspath is needed in case we analyze the "./"
|
|
||||||
folder = os.path.abspath(folder);
|
|
||||||
if basename == 'auto':
|
|
||||||
basename = "_".join(folder.rstrip("/").split("/")[-2:]) + "_"
|
|
||||||
q = data.q if "q" in data else np.arange(data.data.shape[-1])
|
|
||||||
# save one file with all average diffs
|
|
||||||
fname = "%s/%sdiff_av_matrix.txt" % (folder,basename)
|
|
||||||
utils.saveTxt(fname,q,data.data,headerv=data.scan,**kw)
|
|
||||||
# save error bars in the matrix form
|
|
||||||
fname = "%s/%sdiff_av_matrix_err.txt" % (folder,basename)
|
|
||||||
utils.saveTxt(fname,q,data.err,headerv=data.scan,**kw)
|
|
||||||
|
|
||||||
for iscan,scan in enumerate(data.scan):
|
|
||||||
scan = utils.timeToStr(scan) if delayToStr else "%+10.5e" % scan
|
|
||||||
|
|
||||||
# try retreiving info on chi2
|
|
||||||
try:
|
|
||||||
chi2_0 = data.chi2_0[iscan]
|
|
||||||
info_delay = [ "# rep_num : chi2_0 , discarded by chi2masking ?", ]
|
|
||||||
for irep,value in enumerate(chi2_0):
|
|
||||||
info_delay.append( "# %d : %.3f" % (irep,value))
|
|
||||||
if 'chi2' in data.masks: info_delay[-1] += " %s"%str(data.masks['chi2'][iscan][irep])
|
|
||||||
info_delay = "\n".join(info_delay)
|
|
||||||
if info != '': info_delay = "%s\n%s" % (info,info_delay)
|
|
||||||
except AttributeError:
|
|
||||||
info_delay = info
|
|
||||||
|
|
||||||
# save one file per timedelay with average diff (and err)
|
|
||||||
fname = "%s/%sdiff_av_%s.txt" % (folder,basename,scan)
|
|
||||||
if 'mask' in data:
|
|
||||||
tosave = np.vstack( (data.data[iscan],data.err[iscan],
|
|
||||||
data.dataUnmasked[iscan],data.errUnmasked[iscan] ) )
|
|
||||||
columns = 'q diffmask errmask diffnomask errnomask'.split()
|
|
||||||
else:
|
|
||||||
tosave = np.vstack( (data.data[iscan],data.err[iscan] ) )
|
|
||||||
columns = 'q diff err'.split()
|
|
||||||
utils.saveTxt(fname,q,tosave,info=info_delay,columns=columns)
|
|
||||||
|
|
||||||
# save one file per timedelay with all diffs for given delay
|
|
||||||
fname = "%s/%sdiffs_%s.txt" % (folder,basename,scan)
|
|
||||||
utils.saveTxt(fname,q,data.diffsInScanPoint[iscan],info=info_delay,**kw)
|
|
||||||
|
|
||||||
|
|
||||||
def read_diff_av(folder,plot2D=False,save=None):
|
|
||||||
print("Never tested !!!")
|
|
||||||
basename = folder+"/"+"diff_av*"
|
|
||||||
files = glob.glob(basename)
|
|
||||||
files.sort()
|
|
||||||
if len(files) == 0:
|
|
||||||
print("No file found (basename %s)" % basename)
|
|
||||||
return None
|
|
||||||
temp = [os.path.basename(f[:-4]) for f in files]
|
|
||||||
delays = [f.split("_")[-1] for f in temp ]
|
|
||||||
diffav = collections.OrderedDict()
|
|
||||||
diffs = collections.OrderedDict()
|
|
||||||
for d,f in zip(delays,files):
|
|
||||||
data = np.loadtxt(f)
|
|
||||||
diffav[d]=data[:,1]
|
|
||||||
diffs[d] = np.loadtxt(folder+"/diffs_%s.dat"%d)[:,1:]
|
|
||||||
q =data[:,0]
|
|
||||||
t = np.asarray( [mc.strToTime(ti) for ti in delays] )
|
|
||||||
if plot2D:
|
|
||||||
idx = t>0
|
|
||||||
i = np.asarray( diffav.values() )
|
|
||||||
plt.pcolor(np.log10(t[idx]),q,i[idx].T)
|
|
||||||
plt.xlabel(r"$\log_{10}(t)$")
|
|
||||||
plt.ylabel(r"q ($\AA^{-1}$)")
|
|
||||||
it=np.asarray(diffav.values())
|
|
||||||
if save:
|
|
||||||
tosave = np.vstack( (q,it) )
|
|
||||||
header = np.hstack( (len(it),t) )
|
|
||||||
tosave = np.vstack( (header,tosave.T) )
|
|
||||||
np.savetxt(folder + "/all_diffs_av_matrix.txt",tosave)
|
|
||||||
return q,it,diffs,t
|
|
||||||
|
|
@ -1,356 +0,0 @@
|
||||||
from __future__ import print_function,division
|
|
||||||
import os
|
|
||||||
import collections
|
|
||||||
|
|
||||||
# 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
|
|
||||||
|
|
||||||
# use npz files (they can handle more stuff (list of arrays,unicode) than h5py)
|
|
||||||
g_default_ext = '.npz'
|
|
||||||
|
|
||||||
id9.default_extension = g_default_ext
|
|
||||||
g_default_mask = '../masks/fesalen1_run8_careful.edf'
|
|
||||||
g_default_folder = "../fesalen/fesalen1/run%d"
|
|
||||||
|
|
||||||
runsE = collections.OrderedDict()
|
|
||||||
runsT = collections.OrderedDict()
|
|
||||||
runsE[7] = 60; runsT[7] = 135
|
|
||||||
runsE[8] = 60; runsT[8] = 135
|
|
||||||
runsE[9] = 120; runsT[9] = 135
|
|
||||||
runsE[10] = 240; runsT[10] = 135
|
|
||||||
runsE[11] = 30; runsT[11] = 135
|
|
||||||
runsE[12] = 45; runsT[12] = 135
|
|
||||||
runsE[13] = 60; runsT[13] = 120
|
|
||||||
runsE[14] = 120; runsT[14] = 120
|
|
||||||
runsE[15] = 30; runsT[15] = 120
|
|
||||||
|
|
||||||
|
|
||||||
runs135=(11,12,8,9,10)
|
|
||||||
runs120=(15,13,14)
|
|
||||||
|
|
||||||
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 (if needed)
|
|
||||||
if len(logging.getLogger('').handlers)==1:
|
|
||||||
logging.getLogger('').addHandler(console)
|
|
||||||
|
|
||||||
prepareLog()
|
|
||||||
|
|
||||||
def defineCalibrants(dmin=3):
|
|
||||||
c120 = xray.cell.Triclinic(a=10.297,b=10.716,c=13.955,alpha=72.22,beta=70.45,gamma=74.13)
|
|
||||||
p120,_ = c120.reflection_list()
|
|
||||||
c300 = xray.cell.Triclinic(a=10.485,b=11.015,c=14.280,alpha=74.61,beta=69.29,gamma=75.29)
|
|
||||||
p300,_ = c300.reflection_list()
|
|
||||||
return xray.storage.DataStorage( p120=p120, c120=c120,p300=p300,c300=c300)
|
|
||||||
|
|
||||||
def readCalculations(folder="../fesalen/cif",force=False,plot=False):
|
|
||||||
npz = folder + "/fesalen_calc.npz"
|
|
||||||
if not os.path.isfile(npz) or force:
|
|
||||||
ret = dict()
|
|
||||||
for T in ("120LS","200HS","300HS"):
|
|
||||||
# calculated with lambda = 1 Ang
|
|
||||||
theta,i = np.loadtxt( folder+"/fesalen_%s_fwhm0.1.tsv"%T,unpack=True )
|
|
||||||
theta = theta/180*np.pi
|
|
||||||
q = 4*np.pi/1*np.sin(theta/2)
|
|
||||||
key = "T%s"%T
|
|
||||||
ret[ key ] = dict( q = q, i = i )
|
|
||||||
ret = xray.storage.DataStorage(ret)
|
|
||||||
ret.save(folder + "/fesalen_calc.npz" )
|
|
||||||
ret = xray.storage.read(npz)
|
|
||||||
if plot:
|
|
||||||
T = list(ret.keys())
|
|
||||||
T.sort()
|
|
||||||
for i,t in enumerate(T):
|
|
||||||
plt.plot( ret[t].q,ret[t].i, label=t)#,color=plt.cm.Blues((0.3+0.7*i/2)) )
|
|
||||||
plt.grid()
|
|
||||||
plt.legend()
|
|
||||||
plt.xlabel(r"q ($\AA^{-1}$)")
|
|
||||||
return ret
|
|
||||||
|
|
||||||
|
|
||||||
def findCenter(cen=(484,479.5),dist=0.251):
|
|
||||||
files = xray.utils.getEdfFiles("../fesalen/fesalen1/run17/",nFiles=100)
|
|
||||||
img = xray.azav.read(files).mean(axis=0) - 10.
|
|
||||||
wavelength=12.398/18.*1e-10
|
|
||||||
# 1,0,0 peak are the one close to azimuth 0 and are at 0.6596Ang-1 at 120K
|
|
||||||
# and 0.65098 at 300K (run 17 is @ RT)
|
|
||||||
xray.azav.find_center(img,dist=dist,psize=0.0001770834,center=cen,
|
|
||||||
reference=0.65098,wavelength=wavelength)
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
def removeBaseline(folder,qlims=(0.5,2.5),max_iter=30,force=False):
|
|
||||||
if isinstance(folder,int): folder = g_default_folder%folder
|
|
||||||
fname = folder + "/pyfai_1d_nobaseline"+g_default_ext
|
|
||||||
if not os.path.isfile(fname) or force:
|
|
||||||
logging.info("Calling azav with default parameters")
|
|
||||||
data = azav(folder,withBaseline=True)
|
|
||||||
# first peak (around 0.6) is weak and requires special care ...
|
|
||||||
data.q,data.data = xray.utils.removeBackground(data.q,data.data,
|
|
||||||
xlims=qlims,max_iter=max_iter,background_regions=[ (0,0.58) ])
|
|
||||||
data.save(fname)
|
|
||||||
else:
|
|
||||||
data = xray.storage.read(fname)
|
|
||||||
return data
|
|
||||||
|
|
||||||
def azav(folder,nQ=3000,force=False,saveChi=True,mask=g_default_mask,
|
|
||||||
withBaseline=True):
|
|
||||||
if isinstance(folder,int): folder = g_default_folder%folder
|
|
||||||
if withBaseline:
|
|
||||||
data = id9.doFolder_azav(folder,nQ=nQ,force=force,mask=mask,
|
|
||||||
saveChi=saveChi)
|
|
||||||
else:
|
|
||||||
data=removeBaseline(folder)
|
|
||||||
return data
|
|
||||||
|
|
||||||
def datared(folder,monitor=(0.5,4),showPlot=True,errMask=5,chi2Mask=1.5,
|
|
||||||
qlims=(0.5,3),withBaseline=False,forceDatared=False,
|
|
||||||
select=slice(None),**kw):
|
|
||||||
if isinstance(folder,int): folder = g_default_folder%folder
|
|
||||||
data = azav(folder,withBaseline=withBaseline)
|
|
||||||
if withBaseline:
|
|
||||||
diffFile = folder + "/diffs" + g_default_ext
|
|
||||||
else:
|
|
||||||
diffFile = folder + "/diffs_nobaseline" + g_default_ext
|
|
||||||
|
|
||||||
if not os.path.isfile(diffFile) or forceDatared:
|
|
||||||
data,diffs = id9.doFolder_dataRed(data,monitor=monitor,
|
|
||||||
errMask=errMask,chi2Mask=chi2Mask,qlims=qlims,
|
|
||||||
outStorageFile=diffFile)
|
|
||||||
else:
|
|
||||||
diffs = xray.storage.read(diffFile)
|
|
||||||
|
|
||||||
if showPlot:
|
|
||||||
xray.utils.plotdiffs(diffs.q,diffs.data,t=diffs.scan,select=select,
|
|
||||||
absSignal=diffs.dataAbsAvNeg,absSignalScale=10)
|
|
||||||
plt.title(folder + " norm %s" % str(monitor))
|
|
||||||
plt.figure()
|
|
||||||
xray.utils.plotdiffs(diffs.q,diffs.dataAsAbs,t=diffs.scan,
|
|
||||||
select=select)
|
|
||||||
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 _getAmplitudeFromSVD(i,comp=0):
|
|
||||||
u,s,v = np.linalg.svd(i)
|
|
||||||
return u[:,comp]*s[comp]
|
|
||||||
|
|
||||||
def _getAmplitudeFromLinFit(i,ref=1):
|
|
||||||
amp = [mc.linFit(i[ref],ii) for ii in i]
|
|
||||||
return np.asarray(amp)
|
|
||||||
|
|
||||||
def _getAmplitudeFromAbs(i,e):
|
|
||||||
amp = np.abs(i).mean(axis=1)
|
|
||||||
baseline = np.median(e)
|
|
||||||
return amp-baseline
|
|
||||||
|
|
||||||
def _getAmplitudeFromPeak(i):
|
|
||||||
""" get max-min (to be less sensitive to baseline change """
|
|
||||||
amp = np.max(i,axis=1)-np.min(i,axis=1)
|
|
||||||
return amp
|
|
||||||
|
|
||||||
def anaAmplitue(run=12,plot=False,scaleWithE=False,withBaseline=False):
|
|
||||||
folder = g_default_folder % run
|
|
||||||
pars = dict( chi2Mask=1.5,errMask=3,monitor=(0.4,2.5),showPlot=False )
|
|
||||||
data,diffs=datared(folder,**pars,withBaseline=withBaseline)
|
|
||||||
ranges = ( (0.66,0.72), (0.74,0.8), (0.82,0.88) ) #, (0.96,1.02 ) )
|
|
||||||
nPlot = len(ranges)
|
|
||||||
if plot: fig,ax = plt.subplots(nPlot,1,sharex=True)
|
|
||||||
amp = collections.OrderedDict()
|
|
||||||
for num,r in enumerate(ranges):
|
|
||||||
idx = (diffs.q>r[0]) & (diffs.q<r[1])
|
|
||||||
#amplitude = _getAmplitudeFromLinFit( diffs.data[:,idx] )
|
|
||||||
amplitude = _getAmplitudeFromAbs( diffs.data[:,idx],diffs.err[:,idx] )
|
|
||||||
if scaleWithE: amplitude = amplitude/runsE[run]
|
|
||||||
#amplitude = _getAmplitudeFromPeak( diffs.data[:,idx] )
|
|
||||||
amp[r] = amplitude
|
|
||||||
if plot:
|
|
||||||
ax[num].plot(diffs.scan,amplitude,'-o')
|
|
||||||
ax[num].set_title("Range %s"%(str(r)))
|
|
||||||
return diffs.scan,amp
|
|
||||||
|
|
||||||
def _getPeakPos(x,y,e=1):
|
|
||||||
if y.ndim == 1: y = (y,)
|
|
||||||
ret = dict()
|
|
||||||
for yy,ee in zip(y,e):
|
|
||||||
res = xray.peaks.fitPeak(x,yy,err=ee)
|
|
||||||
#plt.figure(); res.plot_fit(); input("ok")
|
|
||||||
pars = res.params.keys()
|
|
||||||
for p in pars:
|
|
||||||
if not p in ret: ret[p] = dict( value = [], err = [] )
|
|
||||||
ret[p]['value'].append( res.params[p].value )
|
|
||||||
ret[p]['err'].append( res.params[p].stderr )
|
|
||||||
for p in pars:
|
|
||||||
ret[p]['value'] = np.asarray(ret[p]['value'])
|
|
||||||
ret[p]['err'] = np.asarray(ret[p]['err'])
|
|
||||||
return xray.storage.DataStorage(ret)
|
|
||||||
|
|
||||||
def anaPeakPos(run=12,plot=False,withBaseline=False):
|
|
||||||
folder = g_default_folder % run
|
|
||||||
pars = dict( chi2Mask=1.5,errMask=3,monitor=(0.4,2.5),showPlot=False )
|
|
||||||
data,diffs=datared(folder,**pars,withBaseline=withBaseline)
|
|
||||||
ranges = ( (0.66,0.72), (0.74,0.8), (0.82,0.88) ) #, (0.96,1.02 ) )
|
|
||||||
nPlot = len(ranges)
|
|
||||||
if plot: fig,ax = plt.subplots(nPlot,1,sharex=True)
|
|
||||||
pars = collections.OrderedDict()
|
|
||||||
for num,r in enumerate(ranges):
|
|
||||||
idx = (diffs.q>r[0]) & (diffs.q<r[1])
|
|
||||||
#amplitude = _getAmplitudeFromLinFit( diffs.data[:,idx] )
|
|
||||||
pars[r] = _getPeakPos( data.q[idx],diffs.dataAsAbs[:,idx],diffs.err[:,idx])
|
|
||||||
if plot:
|
|
||||||
ax[num].axhline( pars[r].center.value[0] )
|
|
||||||
ax[num].plot(diffs.scan,pars[r].center.value,'-o')
|
|
||||||
ax[num].set_title("Range %s"%(str(r)))
|
|
||||||
return diffs.scan,pars
|
|
||||||
|
|
||||||
|
|
||||||
#### PLOTS #####
|
|
||||||
|
|
||||||
def compareT():
|
|
||||||
fig,ax = plt.subplots(2,1,sharex=True)
|
|
||||||
plt.sca( ax[0] )
|
|
||||||
readCalculations(plot=True)
|
|
||||||
plt.sca( ax[1] )
|
|
||||||
t120 = removeBaseline("../fesalen/fesalen1/run16")
|
|
||||||
t300 = removeBaseline("../fesalen/fesalen1/run17")
|
|
||||||
plt.plot(t120.q,t120.data.mean(0),label="Data 120K")
|
|
||||||
plt.plot( [0,.01],[0,.01] ) # just to have a matching color
|
|
||||||
plt.plot(t300.q,t300.data.mean(0),label="Data 300K")
|
|
||||||
|
|
||||||
def showBaseline(N=2):
|
|
||||||
data_bs = azav("../fesalen/fesalen1/run17")
|
|
||||||
data = removeBaseline("../fesalen/fesalen1/run17")
|
|
||||||
idx1 = np.argwhere( data.q[0] == data_bs.q )[0][0]
|
|
||||||
idx2 = np.argwhere( data.q[-1] == data_bs.q )[0][0]
|
|
||||||
bs = data_bs.data[:,idx1:idx2+1]-data.data
|
|
||||||
fig,ax=plt.subplots(2,N,sharex=True,sharey='row',squeeze=False)
|
|
||||||
ax = ax.T
|
|
||||||
for i,a in enumerate(ax):
|
|
||||||
a[0].plot(data_bs.q,data_bs.data[i])
|
|
||||||
a[0].plot(data.q,bs[i])
|
|
||||||
a[1].plot(data.q,data.data[i])
|
|
||||||
|
|
||||||
def showBaselineStability(N=100):
|
|
||||||
# has to be done with cryo on
|
|
||||||
data_bs = azav("../fesalen/fesalen1/run15",withBaseline=True)
|
|
||||||
data = azav("../fesalen/fesalen1/run15",withBaseline=False)
|
|
||||||
fig,ax=plt.subplots(2,2,sharey=True)
|
|
||||||
i = np.arange(data.data.shape[0])
|
|
||||||
ax[0,0].pcolormesh(i,data_bs.q,data_bs.data.T)
|
|
||||||
ax[0,0].set_title("Without baseline subtraction")
|
|
||||||
ax[0,0].set_ylabel(r"q ($\AA^{-1}$)")
|
|
||||||
for x in range(N):
|
|
||||||
ax[0,1].plot(data_bs.data[x],data_bs.q,lw=0.5,color="0.5",alpha=0.1)
|
|
||||||
ax[1,0].pcolormesh(i,data.q,data.data.T)
|
|
||||||
ax[1,0].set_xlabel("Image number")
|
|
||||||
ax[1,0].set_title("With baseline subtraction")
|
|
||||||
ax[1,0].set_ylabel(r"q ($\AA^{-1}$)")
|
|
||||||
for x in range(N):
|
|
||||||
ax[1,1].plot(data.data[x],data.q,lw=0.5,color="0.5",alpha=0.1)
|
|
||||||
ax[0,1].set_xlim(0,400)
|
|
||||||
ax[1,1].set_xlim(0,180)
|
|
||||||
ax[1,1].set_xlabel("Intensity")
|
|
||||||
[a.set_xlim(0,400) for a in ax[:,0]]
|
|
||||||
plt.ylim(0.4,3)
|
|
||||||
|
|
||||||
def comparisonTRbaseline(folder="../fesalen/fesalen1/run9",select=slice(None)):
|
|
||||||
pars = dict( chi2Mask=1.5,errMask=3,monitor=(0.4,2.5),showPlot=False )
|
|
||||||
data,diffs_bs=datared(folder,**pars,withBaseline=True)
|
|
||||||
data,diffs=datared(folder,**pars,withBaseline=False)
|
|
||||||
fig,ax=plt.subplots(2,1,sharex=True)
|
|
||||||
plt.sca(ax[0])
|
|
||||||
xray.utils.plotdiffs(diffs.q,diffs.data,t=diffs.scan,
|
|
||||||
absSignal=diffs.dataAbsAvAll,absSignalScale=10,select=select)
|
|
||||||
plt.title(folder + " norm %s" % str(pars['monitor']))
|
|
||||||
|
|
||||||
plt.sca(ax[1])
|
|
||||||
xray.utils.plotdiffs(diffs_bs.q,diffs_bs.data,t=diffs_bs.scan,
|
|
||||||
absSignal=diffs_bs.dataAbsAvAll,absSignalScale=30,select=select)
|
|
||||||
|
|
||||||
|
|
||||||
def compareAmplitudes(runs=(11,12,8,9,10),scaleWithE=True,withBaseline=False):
|
|
||||||
runs=list(runs)
|
|
||||||
nRuns = len(runs)
|
|
||||||
fig,ax=plt.subplots(nRuns,3,sharex=True,sharey='col',figsize=[12,8])
|
|
||||||
subplots_kw = dict( hspace=0,wspace=0 )
|
|
||||||
plt.subplots_adjust( **subplots_kw )
|
|
||||||
plt.rcParams['font.size']=14
|
|
||||||
for run,arow in zip(runs,ax):
|
|
||||||
t,amp = anaAmplitue(run,plot=False,scaleWithE=scaleWithE,withBaseline=withBaseline)
|
|
||||||
for k,a in zip(amp.keys(),arow):
|
|
||||||
a.axhline(0,ls='--')
|
|
||||||
a.plot(t,amp[k],'o')
|
|
||||||
for r,a in zip(runs,ax[:,0]):
|
|
||||||
label = "run %d (%.0f K)\n%.0f uJ"%(r,runsT[r],runsE[r])
|
|
||||||
a.set_ylabel(label)
|
|
||||||
for interval,a in zip(amp.keys(),ax[0,:]): a.set_title("range %s"%(str(interval)))
|
|
||||||
|
|
||||||
def compareAmplitudes2(runs=(11,12,8,9,10),scaleWithE=False,withBaseline=False):
|
|
||||||
runs=list(runs)
|
|
||||||
nRuns = len(runs)
|
|
||||||
fig,ax=plt.subplots(1,3,sharex=True,figsize=[12,8])
|
|
||||||
subplots_kw = dict( hspace=0,wspace=0 )
|
|
||||||
plt.subplots_adjust( **subplots_kw )
|
|
||||||
plt.rcParams['font.size']=14
|
|
||||||
for a in ax: a.axhline(0,ls='--')
|
|
||||||
for run in runs:
|
|
||||||
t,amp = anaAmplitue(run,plot=False,scaleWithE=scaleWithE,withBaseline=withBaseline)
|
|
||||||
for k,a in zip(amp.keys(),ax):
|
|
||||||
label = "run %d %.0f uJ"%(run,runsE[run])
|
|
||||||
a.plot(t,amp[k],'-o',label=label)
|
|
||||||
# for r,a in zip(runs,ax[:,0]):
|
|
||||||
# a.set_ylabel(label)
|
|
||||||
for interval,a in zip(amp.keys(),ax): a.set_title("range %s"%(str(interval)))
|
|
||||||
|
|
||||||
|
|
||||||
def comparePosition2(runs=(11,12,8,9,10),withBaseline=False,useRatio=True,scaleWithE=False):
|
|
||||||
runs=list(runs)
|
|
||||||
nRuns = len(runs)
|
|
||||||
fig,ax=plt.subplots(1,3,sharex=True,figsize=[12,8],sharey=useRatio)
|
|
||||||
subplots_kw = dict( hspace=0,wspace=0 )
|
|
||||||
plt.subplots_adjust( **subplots_kw )
|
|
||||||
plt.rcParams['font.size']=14
|
|
||||||
if scaleWithE and not ratio:
|
|
||||||
logging.warn("scaleWithE requires useRatio, setting scaleWithE to False")
|
|
||||||
scaleWithE = Fale
|
|
||||||
for run in runs:
|
|
||||||
t,pars = anaPeakPos(run,plot=False,withBaseline=withBaseline)
|
|
||||||
for k,a in zip(pars.keys(),ax):
|
|
||||||
label = "run %d %.0f uJ"%(run,runsE[run])
|
|
||||||
if useRatio:
|
|
||||||
ratio = pars[k].center.value[1:]/pars[k].center.value[0]
|
|
||||||
err = pars[k].center.err[1:]/pars[k].center.value[0]
|
|
||||||
if scaleWithE: ratio = (ratio-1)/runsE[run]+1
|
|
||||||
line = a.errorbar(t[1:],ratio,err,fmt='-o',label=label)[0]
|
|
||||||
else:
|
|
||||||
line = a.plot(t[1:],pars[k].center.value[1:],'-o',label=label)[0]
|
|
||||||
a.axhline(pars[k].center.value[0],color=line.get_color())
|
|
||||||
for interval,a in zip(pars.keys(),ax): a.set_title("range %s"%(str(interval)))
|
|
||||||
if useRatio:
|
|
||||||
ax[0].set_ylabel("Peak Position / Peak Position(t=-5ns)")
|
|
||||||
else:
|
|
||||||
ax[0].set_ylabel("Peak Position ($\AA^{-1}$)")
|
|
||||||
|
|
@ -1,68 +0,0 @@
|
||||||
from __future__ import print_function,division
|
|
||||||
import numpy as np
|
|
||||||
import pylab as plt
|
|
||||||
import os
|
|
||||||
import mcutils as mc
|
|
||||||
import mcutils.xray as xray
|
|
||||||
from mcutils.xray import id9
|
|
||||||
id9 = xray.id9
|
|
||||||
|
|
||||||
g_default_ext = '.npz'
|
|
||||||
# use npz files (they can handle more stuff (list of arrays,unicode) than h5py)
|
|
||||||
id9.default_extension = g_default_ext
|
|
||||||
#id9.default_extension = '.h5'
|
|
||||||
|
|
||||||
def readCalc():
|
|
||||||
fold = "../tiox/calculated_patterns/"
|
|
||||||
names = "alpha beta lam".split()
|
|
||||||
fnames = "alpha500K beta290K lambda".split()
|
|
||||||
calc = dict()
|
|
||||||
for name,fname in zip(names,fnames):
|
|
||||||
q,i=np.loadtxt(fold + "%s.xye.q"%fname,unpack=True)
|
|
||||||
calc[name] = dict( q=q, i=i )
|
|
||||||
return xray.storage.DataStorage(calc)
|
|
||||||
|
|
||||||
calc = readCalc()
|
|
||||||
|
|
||||||
def azav(folder,nQ=1500,force=False,saveChi=True,mask='y>470'):
|
|
||||||
return id9.doFolder_azav(folder,nQ=nQ,force=force,mask=mask,saveChi=saveChi)
|
|
||||||
|
|
||||||
def datared(folder,monitor=(1,5),forceDatared=False,showPlot=True,**kw):
|
|
||||||
azavFile = folder + "/pyfai_1d" + g_default_ext
|
|
||||||
diffFile = folder + "/diffs" + g_default_ext
|
|
||||||
if not os.path.isfile(diffFile) or forceDatared:
|
|
||||||
data,diffs = id9.doFolder_dataRed(folder,monitor=monitor,**kw)
|
|
||||||
else:
|
|
||||||
data = xray.storage.read(azavFile)
|
|
||||||
diffs = xray.storage.read(diffFile)
|
|
||||||
if showPlot:
|
|
||||||
xray.utils.plotdiffs(diffs.q,diffs.data,t=diffs.scan,
|
|
||||||
absSignal=diffs.dataAbsAvAll,absSignalScale=30)
|
|
||||||
plt.plot(calc.lam.q,calc.lam.i/1000+0.1,label='calc')
|
|
||||||
plt.title(folder + " norm %s" % str(monitor))
|
|
||||||
return data,diffs
|
|
||||||
|
|
||||||
def doall(folder,force=False):
|
|
||||||
azav(folder,force=force)
|
|
||||||
return datared(folder,forceDatared=force)
|
|
||||||
|
|
||||||
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)))
|
|
||||||
|
|
||||||
def plotCalc(scale=1):
|
|
||||||
fold = "../tiox/calculated_patterns/"
|
|
||||||
q,i=readtxt(fold + "alpha500K.xye.q")
|
|
||||||
plt.plot(q,i*scale,label="alpha")
|
|
||||||
q,i=readtxt(fold + "beta290K.xye.q")
|
|
||||||
plt.plot(q,i*scale,label="beta")
|
|
||||||
q,i=readtxt(fold + "lambda.xye.q")
|
|
||||||
plt.plot(q,i*scale,label="lambda")
|
|
||||||
|
|
@ -1,85 +0,0 @@
|
||||||
"""
|
|
||||||
module that contains filters and outliers removal procedures
|
|
||||||
most of them return the data array and a dictionary with additional info
|
|
||||||
(parameters, statistics, etc)
|
|
||||||
"""
|
|
||||||
from __future__ import print_function,division
|
|
||||||
from . import utils
|
|
||||||
import logging
|
|
||||||
import statsmodels.robust
|
|
||||||
log = logging.getLogger(__name__) # __name__ is "foo.bar" here
|
|
||||||
|
|
||||||
import numpy as np
|
|
||||||
np.seterr(all='ignore')
|
|
||||||
|
|
||||||
def removeZingers(curves,errs=None,norm='auto',threshold=10,useDerivative=False):
|
|
||||||
""" curves will be normalized internally
|
|
||||||
if errs is None, calculate mad based noise
|
|
||||||
useDerivative for data with trends ..
|
|
||||||
"""
|
|
||||||
# normalize
|
|
||||||
if norm == 'auto':
|
|
||||||
norm = np.nanmean(curves,axis=1)
|
|
||||||
norm = utils.reshapeToBroadcast(norm,curves)
|
|
||||||
|
|
||||||
if useDerivative:
|
|
||||||
data = np.gradient(curves/norn,axis=0)
|
|
||||||
else:
|
|
||||||
data = curves/norm
|
|
||||||
|
|
||||||
median = np.median(data,axis=0)
|
|
||||||
|
|
||||||
# calculate or normalize error
|
|
||||||
if errs is None:
|
|
||||||
errs = statsmodels.robust.mad(data,axis=0)
|
|
||||||
else:
|
|
||||||
errs = errs/norm
|
|
||||||
|
|
||||||
diff = np.abs(data-median)/errs
|
|
||||||
idx = diff > threshold
|
|
||||||
log.debug("Removed %d zingers from %d curves"%(idx.sum(),len(curves)))
|
|
||||||
print("Removed %d zingers from %d curves"%(idx.sum(),len(curves)))
|
|
||||||
if idx.sum()>0:
|
|
||||||
curves[idx]=np.nan
|
|
||||||
#curves = np.ma.MaskedArray(data=curves,mask=idx)
|
|
||||||
return curves
|
|
||||||
|
|
||||||
def filterOutlier(curves,errs=None,norm=None,threshold=10):
|
|
||||||
# normalize
|
|
||||||
if norm == 'auto':
|
|
||||||
norm = np.nanmean(curves,axis=1)
|
|
||||||
norm = utils.reshapeToBroadcast(n,curves)
|
|
||||||
elif norm is None:
|
|
||||||
norm = 1
|
|
||||||
|
|
||||||
curves = curves/norm
|
|
||||||
if errs is None:
|
|
||||||
errs = statsmodels.robust.mad(curves,axis=0)
|
|
||||||
else:
|
|
||||||
errs = errs/norm
|
|
||||||
|
|
||||||
median = np.median(curves)
|
|
||||||
diff = np.abs(curves-median)/errs
|
|
||||||
chi2 = np.sum(diff**2)/len(curves)
|
|
||||||
idx = chi2 < threshold
|
|
||||||
return curves[idx]
|
|
||||||
|
|
||||||
def chi2Filter(diffs,threshold=10):
|
|
||||||
""" Contrary to removeZingers, this removes entire curves """
|
|
||||||
idx_mask = []
|
|
||||||
for iscan in range(len(diffs.diffsInScanPoint)):
|
|
||||||
idx = diffs.chi2_0[iscan] > threshold
|
|
||||||
# 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) )
|
|
||||||
print("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()
|
|
||||||
if "masks_pars" not in data: data['masks_pars'] = dict()
|
|
||||||
data['masks']['chi2'] = idx_mask
|
|
||||||
data['masks_pars']['chi2_threshold'] = threshold
|
|
||||||
return data
|
|
||||||
|
|
||||||
212
xray/id9.py
212
xray/id9.py
|
|
@ -1,212 +0,0 @@
|
||||||
import logging
|
|
||||||
log = logging.getLogger(__name__)
|
|
||||||
|
|
||||||
import os
|
|
||||||
import collections
|
|
||||||
import numpy as np
|
|
||||||
from . import azav
|
|
||||||
from . import dataReduction
|
|
||||||
from . import utils
|
|
||||||
from . import storage
|
|
||||||
from . import filters
|
|
||||||
|
|
||||||
default_extension = ".h5"
|
|
||||||
|
|
||||||
def _conv(x):
|
|
||||||
try:
|
|
||||||
x = float(x)
|
|
||||||
except:
|
|
||||||
x = np.nan
|
|
||||||
return x
|
|
||||||
|
|
||||||
def _readDiagnostic(fname,retry=3):
|
|
||||||
ntry = 0
|
|
||||||
while ntry<retry:
|
|
||||||
try:
|
|
||||||
data = np.genfromtxt(fname,usecols=(2,3),\
|
|
||||||
dtype=None,converters={3: lambda x: _conv(x)},
|
|
||||||
names = ['fname','delay'])
|
|
||||||
return data
|
|
||||||
except Exception as e:
|
|
||||||
log.warn("Could not read diagnostic file, retrying soon,error was %s"%e)
|
|
||||||
ntry += 1
|
|
||||||
# it should not arrive here
|
|
||||||
raise ValueError("Could not read diagnostic file after %d attempts"%retry)
|
|
||||||
|
|
||||||
def readDiagnostic(fname):
|
|
||||||
""" return an ordered dict dictionary of filename; for each key a rounded
|
|
||||||
value of delay is associated """
|
|
||||||
if os.path.isdir(fname): fname += "/diagnostics.log"
|
|
||||||
# try to read diagnostic couple of times
|
|
||||||
data = _readDiagnostic(fname,retry=4)
|
|
||||||
files = data['fname'].astype(str)
|
|
||||||
delays = data['delay']
|
|
||||||
# skip lines that cannot be interpreted as float (like done, etc)
|
|
||||||
idx_ok = np.isfinite( delays )
|
|
||||||
files = files[idx_ok]
|
|
||||||
files = np.asarray( [utils.getBasename(f) for f in files])
|
|
||||||
delays = delays[idx_ok]
|
|
||||||
delays = np.round(delays.astype(float),12)
|
|
||||||
return dict( file = files, scan = delays)
|
|
||||||
|
|
||||||
def _findDark(line):
|
|
||||||
_,value = line.split(":")
|
|
||||||
return float(value)
|
|
||||||
|
|
||||||
def _delayToNum(delay):
|
|
||||||
if delay.decode('ascii') == 'off':
|
|
||||||
delay = -10
|
|
||||||
else:
|
|
||||||
delay=utils.strToTime(delay)
|
|
||||||
return delay
|
|
||||||
|
|
||||||
def findLogFile(folder):
|
|
||||||
files = utils.getFiles(folder,basename='*.log')
|
|
||||||
files.remove(os.path.join(folder,"diagnostics.log"))
|
|
||||||
logfile = files[0]
|
|
||||||
if len(files)>1: log.warn("Found more than one *.log file that is not diagnostics.log: %s"%files)
|
|
||||||
return logfile
|
|
||||||
|
|
||||||
def readLogFile(fnameOrFolder,subtractDark=False,skip_first=0,asDataStorage=True,last=None):
|
|
||||||
""" read id9 style logfile """
|
|
||||||
if os.path.isdir(fnameOrFolder):
|
|
||||||
fname = findLogFile(fnameOrFolder)
|
|
||||||
else:
|
|
||||||
fname = fnameOrFolder
|
|
||||||
f = open(fname,"r")
|
|
||||||
lines = f.readlines()
|
|
||||||
f.close()
|
|
||||||
lines = [line.strip() for line in lines]
|
|
||||||
darks = {}
|
|
||||||
for line in lines:
|
|
||||||
if line.find("pd1 dark/sec")>=0: darks['pd1ic'] = _findDark(line)
|
|
||||||
if line.find("pd2 dark/sec")>=0: darks['pd2ic'] = _findDark(line)
|
|
||||||
if line.find("pd3 dark/sec")>=0: darks['pd3ic'] = _findDark(line)
|
|
||||||
if line.find("pd4 dark/sec")>=0: darks['pd4ic'] = _findDark(line)
|
|
||||||
for iline,line in enumerate(lines):
|
|
||||||
if line.lstrip()[0] != "#": break
|
|
||||||
data=np.genfromtxt(fname,skip_header=iline-1,names=True,comments="%",dtype=None,converters = {'delay': lambda s: _delayToNum(s)})
|
|
||||||
if subtractDark:
|
|
||||||
for diode in ['pd1ic','pd2ic','pd3ic','pd4ic']:
|
|
||||||
if diode in darks: data[diode]=data[diode]-darks[diode]*data['timeic']
|
|
||||||
data = data[skip_first:last]
|
|
||||||
if asDataStorage:
|
|
||||||
# rstrip("_") is used to clean up last _ that appera for some reason in file_
|
|
||||||
data = storage.DataStorage( dict((name.rstrip("_"),data[name]) for name in data.dtype.names ) )
|
|
||||||
data.file = data.file.astype(str)
|
|
||||||
return data
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
def doFolder_azav(folder,nQ=1500,files='*.edf*',force=False,mask=None,
|
|
||||||
saveChi=True,poni='pyfai.poni',storageFile='auto',dark=9.9,dezinger=None,
|
|
||||||
qlims=(0,10),removeBack=False,removeBack_kw=dict(),skip_first=0,
|
|
||||||
last=None):
|
|
||||||
""" very small wrapper around azav.doFolder, essentially just reading
|
|
||||||
the id9 logfile or diagnostics.log """
|
|
||||||
|
|
||||||
try:
|
|
||||||
loginfo = readLogFile(folder,skip_first=skip_first,last=last)
|
|
||||||
except Exception as e:
|
|
||||||
log.warn("Could not read log file, trying to read diagnostics.log")
|
|
||||||
loginfo = readDiagnostic(folder)
|
|
||||||
if storageFile == 'auto' : storageFile = folder + "/" + "pyfai_1d" + default_extension
|
|
||||||
|
|
||||||
data = azav.doFolder(folder,files=files,nQ=nQ,force=force,mask=mask,
|
|
||||||
saveChi=saveChi,poni=poni,storageFile=storageFile,logDict=loginfo,dark=dark,save=False,dezinger=dezinger)
|
|
||||||
#try:
|
|
||||||
# if removeBack is not None:
|
|
||||||
# _,data.data = azav.removeBackground(data,qlims=qlims,**removeBack_kw)
|
|
||||||
#except Exception as e:
|
|
||||||
# log.error("Could not remove background, error was %s"%(str(e)))
|
|
||||||
|
|
||||||
# idx = utils.findSlice(data.q,qlims)
|
|
||||||
# n = np.nanmean(data.data[:,idx],axis=1)
|
|
||||||
# data.norm_range = qlims
|
|
||||||
# data.norm = n
|
|
||||||
# n = utils.reshapeToBroadcast(n,data.data)
|
|
||||||
# data.data_norm = data.data/n
|
|
||||||
|
|
||||||
data.save(storageFile)
|
|
||||||
|
|
||||||
|
|
||||||
return data
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
def doFolder_dataRed(azavStorage,monitor=None,funcForAveraging=np.nanmean,
|
|
||||||
qlims=None,outStorageFile='auto',reference='min'):
|
|
||||||
""" azavStorage if a DataStorage instance or the filename to read
|
|
||||||
monitor : normalization vector that can be given as
|
|
||||||
1. numpy array
|
|
||||||
2. a list (interpreted as q-range normalization)
|
|
||||||
3. a string to look for as key in the log, e.g.
|
|
||||||
monitor="pd2ic" would reult in using
|
|
||||||
azavStorage.log.pd2ic
|
|
||||||
|
|
||||||
"""
|
|
||||||
|
|
||||||
if isinstance(azavStorage,storage.DataStorage):
|
|
||||||
data = azavStorage
|
|
||||||
folder = azavStorage.folder
|
|
||||||
elif os.path.isfile(azavStorage):
|
|
||||||
folder = os.path.dirname(azavStorage)
|
|
||||||
data = storage.DataStorage(azavStorage)
|
|
||||||
else:
|
|
||||||
# assume is just a folder name
|
|
||||||
folder = azavStorage
|
|
||||||
azavStorage = folder + "/pyfai_1d" + default_extension
|
|
||||||
data = storage.DataStorage(azavStorage)
|
|
||||||
|
|
||||||
#assert data.q.shape[0] == data.data.shape[1] == data.err.shape[1]
|
|
||||||
if qlims is not None:
|
|
||||||
idx = (data.q>qlims[0]) & (data.q<qlims[1])
|
|
||||||
data.data = data.data[:,idx]
|
|
||||||
data.err = data.err[:,idx]
|
|
||||||
data.q = data.q[idx]
|
|
||||||
|
|
||||||
if isinstance(monitor,str): monitor = data['log'][monitor]
|
|
||||||
|
|
||||||
# calculate differences
|
|
||||||
diffs = dataReduction.calcTimeResolvedSignal(data.log.delay,data.data,
|
|
||||||
err=data.err,q=data.q,reference=reference,monitor=monitor,
|
|
||||||
funcForAveraging=funcForAveraging)
|
|
||||||
|
|
||||||
# save txt and npz file
|
|
||||||
dataReduction.saveTxt(folder,diffs,info=data.pyfai_info)
|
|
||||||
if outStorageFile == 'auto':
|
|
||||||
outStorageFile = folder + "/diffs" + default_extension
|
|
||||||
diffs.save(outStorageFile)
|
|
||||||
|
|
||||||
return data,diffs
|
|
||||||
|
|
||||||
def doFolder(folder,azav_kw = dict(), datared_kw = dict(),online=True, retryMax=20,force=False):
|
|
||||||
import matplotlib.pyplot as plt
|
|
||||||
if folder == "./": folder = os.path.abspath(folder)
|
|
||||||
fig = plt.figure()
|
|
||||||
lastNum = None
|
|
||||||
keepGoing = True
|
|
||||||
lines = None
|
|
||||||
retryNum = 0
|
|
||||||
if online: print("Press Ctrl+C to stop")
|
|
||||||
while keepGoing and retryNum < retryMax:
|
|
||||||
try:
|
|
||||||
data = doFolder_azav(folder,**azav_kw)
|
|
||||||
# check if there are new data
|
|
||||||
if lastNum is None or lastNum<data.data.shape[0]:
|
|
||||||
data,diffs = doFolder_dataRed(data,**datared_kw)
|
|
||||||
if lines is None or len(lines) != diffs.data.shape[0]:
|
|
||||||
lines,_ = utils.plotdiffs(diffs,fig=fig,title=folder)
|
|
||||||
else:
|
|
||||||
utils.updateLines(lines,diffs.data)
|
|
||||||
plt.draw()
|
|
||||||
lastNum = data.data.shape[0]
|
|
||||||
retryNum = 0
|
|
||||||
else:
|
|
||||||
retryNum += 1
|
|
||||||
plt.pause(30)
|
|
||||||
if force: force=False; # does not make sense to have always True ...
|
|
||||||
except KeyboardInterrupt:
|
|
||||||
keepGoing = False
|
|
||||||
if not online: keepGoing = False
|
|
||||||
return data,diffs
|
|
||||||
210
xray/mask.py
210
xray/mask.py
|
|
@ -1,210 +0,0 @@
|
||||||
from __future__ import print_function
|
|
||||||
import sys
|
|
||||||
if sys.version_info.major == 2: input=raw_input
|
|
||||||
import logging
|
|
||||||
log = logging.getLogger(__name__)
|
|
||||||
|
|
||||||
import os
|
|
||||||
import collections
|
|
||||||
import numpy as np
|
|
||||||
import matplotlib.pyplot as plt
|
|
||||||
from matplotlib.path import Path
|
|
||||||
|
|
||||||
maskComponent = collections.namedtuple('maskComponent',['operation','geometry','vertices'])
|
|
||||||
|
|
||||||
def _rectangleToMask(X,Y,vertices):
|
|
||||||
( (x1,y1), (x2,y2) ) = vertices
|
|
||||||
if x1>x2: x1,x2=x2,x1
|
|
||||||
if y1>y2: y1,y2=y2,y1
|
|
||||||
return (X>x1) & (X<x2) & ( Y>y1) & (Y<y2)
|
|
||||||
|
|
||||||
def _circleToMask(X,Y,vertices):
|
|
||||||
c,p = vertices
|
|
||||||
r = np.sqrt( (p[0]-c[0])**2 + (p[1]-c[1])**2 )
|
|
||||||
d = np.sqrt((X-c[0])**2+(Y-c[1])**2)
|
|
||||||
return d<r
|
|
||||||
|
|
||||||
def _polygonToMask(X,Y,vertices):
|
|
||||||
points = np.vstack((X.flatten(),Y.flatten())).T
|
|
||||||
path = Path(vertices)
|
|
||||||
grid = path.contains_points(points)
|
|
||||||
return grid.reshape(X.shape)
|
|
||||||
|
|
||||||
class MyMask(object):
|
|
||||||
def __init__(self,img=None):
|
|
||||||
self.comp = []
|
|
||||||
self.img = img
|
|
||||||
self.mask = None
|
|
||||||
self._cache = None
|
|
||||||
|
|
||||||
def _define_component(self,operation,geometry,*vertices):
|
|
||||||
#print("define_comp",type(vertices),vertices)
|
|
||||||
if geometry == 'circle' and len(vertices) == 3:
|
|
||||||
xcen,ycen,radius = vertices
|
|
||||||
vertices = ( (xcen,ycen), (xcen+radius,ycen) )
|
|
||||||
if geometry == 'rectangle' and len(vertices) == 4:
|
|
||||||
vertices = ( (vertices[0],vertices[1]),(vertices[2],vertices[3]) )
|
|
||||||
# make sure vertices tuples
|
|
||||||
if isinstance(vertices,list):
|
|
||||||
vertices = [ (v[0],v[1]) for v in vertices ]
|
|
||||||
vertices = tuple(vertices)
|
|
||||||
a = dict( vertices = None )
|
|
||||||
self.comp.append( maskComponent(operation=operation,vertices=vertices,geometry=geometry) )
|
|
||||||
|
|
||||||
|
|
||||||
def addCircle(self,*vertices): self._define_component( 'add', 'circle', *vertices )
|
|
||||||
def subtractCircle(self,*vertices): self._define_component( 'subtract', 'circle', *vertices )
|
|
||||||
|
|
||||||
def addRectangle(self,*vertices): self._define_component( 'add','rectangle', *vertices)
|
|
||||||
def subtractRectangle(self,*vertices): self._define_component( 'subtract','rectangle',*vertices)
|
|
||||||
|
|
||||||
|
|
||||||
def addPolygon(self,*vertices): self._define_component( 'add','polygon',*vertices)
|
|
||||||
def subtractPolygon(self,*vertices): self._define_component( 'subtract','polygon',*vertices)
|
|
||||||
|
|
||||||
|
|
||||||
def getMask(self,shape=None):
|
|
||||||
if shape is None and self.img is not None: shape = self.img.shape
|
|
||||||
if shape is None and self.img is None: shape = self._cache['shape']
|
|
||||||
|
|
||||||
if self._cache is None: self._cache = dict( shape = shape )
|
|
||||||
# reset cache if shape does not match
|
|
||||||
if shape != self._cache['shape']:
|
|
||||||
self._cache = dict( shape = shape )
|
|
||||||
X,Y = np.meshgrid ( range(shape[1]),range(shape[0]) )
|
|
||||||
for component in self.comp:
|
|
||||||
if component not in self._cache:
|
|
||||||
if component.geometry == 'circle':
|
|
||||||
mask = _circleToMask( X,Y,component.vertices )
|
|
||||||
elif component.geometry == 'rectangle':
|
|
||||||
mask = _rectangleToMask( X,Y,component.vertices )
|
|
||||||
elif component.geometry == 'polygon':
|
|
||||||
mask = _polygonToMask( X,Y,component.vertices )
|
|
||||||
else:
|
|
||||||
raise ValueError("Mask type %s not recongnized"%component.geometry)
|
|
||||||
self._cache[component] = mask
|
|
||||||
mask = np.zeros(shape,dtype=np.bool)
|
|
||||||
for comp in self.comp:
|
|
||||||
m = self._cache[ comp ]
|
|
||||||
if (comp.operation == "add"):
|
|
||||||
mask[m] = True
|
|
||||||
else:
|
|
||||||
mask[m] = False
|
|
||||||
self.mask = mask
|
|
||||||
return mask
|
|
||||||
|
|
||||||
def getMatplotlibMask(self,shape=None):
|
|
||||||
mask = self.getMask(shape=shape)
|
|
||||||
# convert
|
|
||||||
mpl_mask = np.zeros( (mask.shape[0],mask.shape[1],4) )
|
|
||||||
mpl_mask[:,:,:3] = 0.5; # gray color
|
|
||||||
mpl_mask[:,:,3] = mask/2; # give some transparency
|
|
||||||
return mpl_mask
|
|
||||||
|
|
||||||
def save(self,fname,inverted=False):
|
|
||||||
import fabio
|
|
||||||
mask = self.mask
|
|
||||||
if (inverted): mask = ~mask
|
|
||||||
i=fabio.edfimage.edfimage(mask.astype(np.uint8)); # edf does not support bool
|
|
||||||
i.save(fname)
|
|
||||||
|
|
||||||
def snap(point,shape,snapRange=20):
|
|
||||||
snapped = list(point)
|
|
||||||
if snapped[0] < snapRange: snapped[0] = 0
|
|
||||||
if snapped[0] > shape[1]-snapRange: snapped[0] = shape[1]
|
|
||||||
if snapped[1] < snapRange: snapped[1] = 0
|
|
||||||
if snapped[1] > shape[0]-snapRange: snapped[1] = shape[0]
|
|
||||||
return tuple(snapped)
|
|
||||||
|
|
||||||
def getPoints(N=1,shape=(100,100),snapRange=0):
|
|
||||||
if N<1: print('Right click cancels last point, middle click ends the polygon')
|
|
||||||
c = plt.ginput(N)
|
|
||||||
c = [ snap(point,shape,snapRange=snapRange) for point in c ]
|
|
||||||
if len(c) == 1: c = c[0]
|
|
||||||
return c
|
|
||||||
|
|
||||||
def makeMaskGui(img,snapRange=60):
|
|
||||||
""" snapRange controls border snapping (in pixels, use <= 0 to disable """
|
|
||||||
mask = MyMask(img)
|
|
||||||
ans='ok'
|
|
||||||
while (ans != 'done'):
|
|
||||||
plt.imshow(img)
|
|
||||||
plt.clim(np.percentile(img,(2,98)))
|
|
||||||
plt.imshow(mask.getMatplotlibMask())
|
|
||||||
plt.pause(0.01)
|
|
||||||
ans = input("What's next p/P/c/C/r/R/done? (capitals = subtract)")
|
|
||||||
if ans == "c":
|
|
||||||
print("Adding circle, click on center then another point to define radius")
|
|
||||||
vertices = getPoints(N=2,shape=img.shape,snapRange=snapRange)
|
|
||||||
mask.addCircle(*vertices)
|
|
||||||
if ans == "C":
|
|
||||||
print("Subtracting circle, click on center then another point to define radius")
|
|
||||||
vertices = getPoints(N=2,shape=img.shape,snapRange=snapRange)
|
|
||||||
mask.subtractCircle(*vertices)
|
|
||||||
if ans == "r":
|
|
||||||
print("Adding rectangle, click on one corner and then on the opposite one")
|
|
||||||
vertices = getPoints(N=2,shape=img.shape,snapRange=snapRange)
|
|
||||||
mask.addRectangle(*vertices)
|
|
||||||
if ans == "R":
|
|
||||||
print("Subtracting rectangle, click on one corner and then on the opposite one")
|
|
||||||
vertices = getPoints(N=2,shape=img.shape,snapRange=snapRange)
|
|
||||||
mask.subtractRectangle(*vertices)
|
|
||||||
if ans == 'p':
|
|
||||||
print("Adding polygon")
|
|
||||||
vertices = getPoints(N=-1,shape=img.shape,snapRange=snapRange)
|
|
||||||
mask.addPolygon(*vertices)
|
|
||||||
if ans == 'P':
|
|
||||||
print("Subtracting polygon")
|
|
||||||
vertices = getPoints(N=-1,shape=img.shape,snapRange=snapRange)
|
|
||||||
mask.subtractPolygon(*vertices)
|
|
||||||
|
|
||||||
plt.imshow(mask.getMatplotlibMask())
|
|
||||||
plt.pause(0.01)
|
|
||||||
fname = input("Enter a valid filename (ext .edf or .npy) if you want to save the mask (empty otherwise)")
|
|
||||||
try:
|
|
||||||
if fname != '':
|
|
||||||
ext = os.path.splitext(fname)[1]
|
|
||||||
if ext == '.edf':
|
|
||||||
mask.save(fname)
|
|
||||||
elif ext == '.npy':
|
|
||||||
np.save(fname,mask.getMask())
|
|
||||||
except Exception as e:
|
|
||||||
log.error("Error in saving mask")
|
|
||||||
log.error(e)
|
|
||||||
finally:
|
|
||||||
return mask
|
|
||||||
|
|
||||||
def maskBorder(width,shape):
|
|
||||||
mask = np.zeros(shape,dtype=bool)
|
|
||||||
mask[ :width , : ] = True
|
|
||||||
mask[ -width: , : ] = True
|
|
||||||
mask[ : , :width ] = True
|
|
||||||
mask[ : , -width: ] = True
|
|
||||||
return mask
|
|
||||||
|
|
||||||
def maskCenterLines(width,shape):
|
|
||||||
mask = np.zeros(shape,dtype=bool)
|
|
||||||
if isinstance(width,int): width = (width,width)
|
|
||||||
c0 = int(shape[0]/2)
|
|
||||||
c1 = int(shape[1]/2)
|
|
||||||
w0 = int(width[0]/2)
|
|
||||||
w1 = int(width[1]/2)
|
|
||||||
mask[ c0-w0:c0+w0 , : ] = True
|
|
||||||
mask[ : , c1-w1:c1+w1 ] = True
|
|
||||||
return mask
|
|
||||||
|
|
||||||
|
|
||||||
def test(shape=(1000,2000)):
|
|
||||||
mask = MyMask()
|
|
||||||
mask.addCircle(400,300,250)
|
|
||||||
mask.subtractCircle(400,300,150)
|
|
||||||
mask.addRectangle(350,250,1500,700)
|
|
||||||
plt.imshow( mask.getMask(shape) )
|
|
||||||
return mask
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
if __name__ == "__main__":
|
|
||||||
test()
|
|
||||||
plt.show()
|
|
||||||
ans=input("Enter to finish")
|
|
||||||
|
|
@ -1,32 +0,0 @@
|
||||||
import lmfit
|
|
||||||
import numpy as np
|
|
||||||
import logging as log
|
|
||||||
|
|
||||||
pv = lmfit.models.PseudoVoigtModel()
|
|
||||||
|
|
||||||
def fitPeak(x,y,err=1,autorange=False):
|
|
||||||
if isinstance(err,np.ndarray):
|
|
||||||
if np.all(err==0):
|
|
||||||
err = 1
|
|
||||||
log.warn("Asked to fit peak but all errors are zero, forcing them to 1")
|
|
||||||
elif np.isfinite(err).sum() != len(err):
|
|
||||||
idx = np.isfinite(err)
|
|
||||||
x = x[idx]
|
|
||||||
y = y[idx]
|
|
||||||
err = err[idx]
|
|
||||||
log.warn("Asked to fit peak but some errorbars are infinite or nans,\
|
|
||||||
excluding those points")
|
|
||||||
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,weights=1/err,params=pars)
|
|
||||||
return ret
|
|
||||||
346
xray/utils.py
346
xray/utils.py
|
|
@ -1,346 +0,0 @@
|
||||||
from __future__ import print_function,division
|
|
||||||
|
|
||||||
import logging
|
|
||||||
log = logging.getLogger(__name__) # __name__ is "foo.bar" here
|
|
||||||
|
|
||||||
import numpy as np
|
|
||||||
np.seterr(all='ignore')
|
|
||||||
import os
|
|
||||||
import glob
|
|
||||||
import pathlib
|
|
||||||
import re
|
|
||||||
import numbers
|
|
||||||
from . import storage as storage
|
|
||||||
|
|
||||||
try:
|
|
||||||
import matplotlib.pyplot as plt
|
|
||||||
except ImportError:
|
|
||||||
log.warn("Can't import matplotlib !")
|
|
||||||
|
|
||||||
_time_regex = re.compile( "(-?\d+\.?\d*(?:ps|ns|us|ms)?)")
|
|
||||||
_timeInStr_regex = re.compile("_(-?\d+\.?\d*(?:ps|ns|us|ms)?)")
|
|
||||||
|
|
||||||
def getFiles(folder,basename="*.edf*",nFiles=None):
|
|
||||||
files = glob.glob(folder + "/" + basename)
|
|
||||||
files.sort()
|
|
||||||
if nFiles is not None: files = files[:nFiles]
|
|
||||||
return files
|
|
||||||
|
|
||||||
def getEdfFiles(folder,**kw):
|
|
||||||
return getFiles(folder,basename="*.edf*",**kw)
|
|
||||||
|
|
||||||
def getDelayFromString(string) :
|
|
||||||
match = _timeInStr_regex.search(string)
|
|
||||||
return match and match.group(1) or None
|
|
||||||
|
|
||||||
_time_regex = re.compile("(-?\d+\.?\d*)((?:s|fs|ms|ns|ps|us)?)")
|
|
||||||
|
|
||||||
def strToTime(delay) :
|
|
||||||
if isinstance(delay,bytes): delay = delay.decode('ascii')
|
|
||||||
_time2value = dict( fs = 1e-15, ps = 1e-12, ns = 1e-9, us = 1e-6, ms = 1e-3, s = 1)
|
|
||||||
|
|
||||||
match = _time_regex.search(delay)
|
|
||||||
if match:
|
|
||||||
n,t = float(match.group(1)),match.group(2)
|
|
||||||
value = _time2value.get(t,1)
|
|
||||||
return n*value
|
|
||||||
else:
|
|
||||||
return None
|
|
||||||
|
|
||||||
def timeToStr(delay,fmt="%+.0f"):
|
|
||||||
a_delay = abs(delay)
|
|
||||||
if a_delay >= 1:
|
|
||||||
ret = fmt % delay + "s"
|
|
||||||
elif 1e-3 <= a_delay < 1:
|
|
||||||
ret = fmt % (delay*1e3) + "ms"
|
|
||||||
elif 1e-6 <= a_delay < 1e-3:
|
|
||||||
ret = fmt % (delay*1e6) + "us"
|
|
||||||
elif 1e-9 <= a_delay < 1e-6:
|
|
||||||
ret = fmt % (delay*1e9) + "ns"
|
|
||||||
elif 1e-12 <= a_delay < 1e-9:
|
|
||||||
ret = fmt % (delay*1e12) + "ps"
|
|
||||||
elif 1e-15 <= a_delay < 1e-12:
|
|
||||||
ret = fmt % (delay*1e12) + "fs"
|
|
||||||
elif 1e-18 <= a_delay < 1e-15:
|
|
||||||
ret = fmt % (delay*1e12) + "as"
|
|
||||||
else:
|
|
||||||
ret = str(delay) +"s"
|
|
||||||
return ret
|
|
||||||
|
|
||||||
def removeExt(fname):
|
|
||||||
""" special remove extension meant to work with compressed files.edf and .edf.gz files """
|
|
||||||
if fname[-3:] == ".gz": fname = fname[:-3]
|
|
||||||
return os.path.splitext(fname)[0]
|
|
||||||
|
|
||||||
def getBasename(fname):
|
|
||||||
return os.path.basename(removeExt(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(*args,x=None,plot=True,showTrend=True,title=None,clim='auto',fig=None):
|
|
||||||
if isinstance(args[0],storage.DataStorage):
|
|
||||||
q = args[0].q; data=args[0].data;
|
|
||||||
if title is None: title = args[0].folder
|
|
||||||
else:
|
|
||||||
q,data = args
|
|
||||||
if not (plot or showTrend): return
|
|
||||||
if x is None: x = np.arange(data.shape[0])
|
|
||||||
if clim == 'auto': clim = np.nanpercentile(data,(1.5,98.5))
|
|
||||||
one_plot = showTrend or plot
|
|
||||||
two_plot = showTrend and plot
|
|
||||||
if one_plot and not two_plot:
|
|
||||||
if fig is None:
|
|
||||||
fig,ax = plt.subplots(1,1)
|
|
||||||
else:
|
|
||||||
fig.clear()
|
|
||||||
ax = fig.axes
|
|
||||||
if two_plot:
|
|
||||||
if fig is None:
|
|
||||||
fig,ax = plt.subplots(2,1,sharex=True)
|
|
||||||
else:
|
|
||||||
ax = fig.axes
|
|
||||||
|
|
||||||
ax = np.atleast_1d(ax)
|
|
||||||
if showTrend:
|
|
||||||
plt.sca(ax[1])
|
|
||||||
plt.pcolormesh(q,x,data)
|
|
||||||
plt.ylabel("image number, 0 being older")
|
|
||||||
plt.xlabel(r"q ($\AA^{-1}$)")
|
|
||||||
plt.clim( *clim )
|
|
||||||
if plot:
|
|
||||||
ax[0].plot(q,np.nanmean(data,axis=0))
|
|
||||||
if (plot or showTrend) and title is not None:
|
|
||||||
plt.title(title)
|
|
||||||
|
|
||||||
def volumeFraction(concentration=1,molWeight=17,density=1.347):
|
|
||||||
""" molWeight is kDa
|
|
||||||
concentration in mM
|
|
||||||
density g/ml
|
|
||||||
"""
|
|
||||||
concentration_mg_ml = concentration*molWeight
|
|
||||||
volume_fraction = concentration_mg_ml/density/1e3
|
|
||||||
return volume_fraction
|
|
||||||
|
|
||||||
def plotdiffs(*args,select=None,err=None,absSignal=None,absSignalScale=10,
|
|
||||||
showErr=False,cmap=plt.cm.jet,fig=None,title=None):
|
|
||||||
# this selection trick done in this way allows to keep the same colors when
|
|
||||||
# subselecting (because I do not change the size of diffs)
|
|
||||||
if isinstance(args[0],storage.DataStorage):
|
|
||||||
q = args[0].q; t = args[0].scan; err = args[0].err
|
|
||||||
diffs = args[0].data
|
|
||||||
diffs_abs = args[0].dataAbsAvScanPoint
|
|
||||||
else:
|
|
||||||
q,diffs,t = args
|
|
||||||
diffs_abs = None
|
|
||||||
if select is not None:
|
|
||||||
indices = range(*select.indices(t.shape[0]))
|
|
||||||
else:
|
|
||||||
indices = range(len(t))
|
|
||||||
|
|
||||||
if fig is None: fig = plt.gcf()
|
|
||||||
# fig.clear()
|
|
||||||
|
|
||||||
lines_diff = []
|
|
||||||
lines_abs = []
|
|
||||||
if absSignal is not None:
|
|
||||||
line = plt.plot(q,absSignal/absSignalScale,lw=3,
|
|
||||||
color='k',label="absSignal/%s"%str(absSignalScale))[0]
|
|
||||||
lines.append(line)
|
|
||||||
for linenum,idiff in enumerate(indices):
|
|
||||||
color = cmap(idiff/(len(diffs)-1))
|
|
||||||
label = timeToStr(t[idiff])
|
|
||||||
kw = dict( color = color, label = label )
|
|
||||||
if err is not None and showErr:
|
|
||||||
line = plt.errorbar(q,diffs[idiff],err[idiff],**kw)[0]
|
|
||||||
lines_diff.append(line)
|
|
||||||
else:
|
|
||||||
line = plt.plot(q,diffs[idiff],**kw)[0]
|
|
||||||
lines_diff.append(line)
|
|
||||||
if diffs_abs is not None:
|
|
||||||
line = plt.plot(q,diffs_abs[idiff],color=color)[0]
|
|
||||||
lines_abs.append(line)
|
|
||||||
if title is not None: fig.axes[0].set_title(title)
|
|
||||||
legend = plt.legend(loc=4)
|
|
||||||
plt.grid()
|
|
||||||
plt.xlabel(r"q ($\AA^{-1}$)")
|
|
||||||
# we will set up a dict mapping legend line to orig line, and enable
|
|
||||||
# picking on the legend line
|
|
||||||
lined = dict()
|
|
||||||
for legline, origline in zip(legend.get_lines(), lines_diff):
|
|
||||||
legline.set_picker(5) # 5 pts tolerance
|
|
||||||
lined[legline] = origline
|
|
||||||
|
|
||||||
def onpick(event):
|
|
||||||
# on the pick event, find the orig line corresponding to the
|
|
||||||
# legend proxy line, and toggle the visibility
|
|
||||||
legline = event.artist
|
|
||||||
origline = lined[legline]
|
|
||||||
vis = not origline.get_visible()
|
|
||||||
origline.set_visible(vis)
|
|
||||||
# Change the alpha on the line in the legend so we can see what lines
|
|
||||||
# have been toggled
|
|
||||||
if vis:
|
|
||||||
legline.set_alpha(1.0)
|
|
||||||
else:
|
|
||||||
legline.set_alpha(0.2)
|
|
||||||
fig = plt.gcf()
|
|
||||||
fig.canvas.draw()
|
|
||||||
|
|
||||||
fig.canvas.mpl_connect('pick_event', onpick)
|
|
||||||
return lines_diff,lines_abs
|
|
||||||
|
|
||||||
def updateLines(lines,data):
|
|
||||||
for l,d in zip(lines,data):
|
|
||||||
l.set_ydata(d)
|
|
||||||
|
|
||||||
|
|
||||||
#def getScan
|
|
||||||
|
|
||||||
def saveTxt(fname,q,data,headerv=None,info=None,overwrite=True,columns=''):
|
|
||||||
""" Write data to file 'fname' in text format.
|
|
||||||
Inputs:
|
|
||||||
q = x vector
|
|
||||||
data = one or 2D array (first axis is q)
|
|
||||||
info = dictionary (saved as '# key : value') or string
|
|
||||||
headerv = vector to be used as header or string
|
|
||||||
"""
|
|
||||||
if os.path.isfile(fname) and not overwrite:
|
|
||||||
log.warn("File %s exists, returning",fname)
|
|
||||||
return
|
|
||||||
if isinstance(info,dict):
|
|
||||||
header = [ "# %s : %s" %(k,str(v)) for (k,v) in info.items() ]
|
|
||||||
header = "\n".join(header); # skip first #, will be added by np
|
|
||||||
elif isinstance(info,str):
|
|
||||||
header = info
|
|
||||||
else:
|
|
||||||
header = ""
|
|
||||||
if isinstance(headerv,str): header += "\n%s" % headerv
|
|
||||||
if data.ndim == 1:
|
|
||||||
x = np.vstack( (q,data) )
|
|
||||||
elif data.ndim == 2:
|
|
||||||
x = np.vstack( (q,data) )
|
|
||||||
if headerv is not None:
|
|
||||||
headerv = np.concatenate(( (data.shape[1],),headerv))
|
|
||||||
x = np.hstack( (headerv[:,np.newaxis],x) )
|
|
||||||
if columns != '':
|
|
||||||
s = "#" + " ".join( [str(c).center(12) for c in columns] )
|
|
||||||
header = header + "\n" + s if header != '' else s
|
|
||||||
np.savetxt(fname,x.T,fmt="%+10.5e",header=header,comments='')
|
|
||||||
|
|
||||||
def reshapeToBroadcast(what,ref):
|
|
||||||
""" expand the 1d array 'what' to allow broadbasting to match
|
|
||||||
multidimentional array 'ref'. The two arrays have to same the same
|
|
||||||
dimensions along the first axis
|
|
||||||
"""
|
|
||||||
if what.shape == ref.shape: return what
|
|
||||||
assert what.shape[0] == ref.shape[0]
|
|
||||||
shape = [ref.shape[0],] + [1,]*(ref.ndim-1)
|
|
||||||
return what.reshape(shape)
|
|
||||||
|
|
||||||
def radToQ(theta,**kw):
|
|
||||||
""" theta is the scattering angle (theta_out - theta_in);
|
|
||||||
kw should be have either E or wavelength
|
|
||||||
it returns the scattering vector in the units of wavelength """
|
|
||||||
# Energy or wavelength should be in kw
|
|
||||||
assert "E" in kw or "wavelength" in kw
|
|
||||||
# but not both
|
|
||||||
assert not ("E" in kw and "wavelength" in kw)
|
|
||||||
if "E" in kw: kw["wavelength"] = 12.398/kw["E"]
|
|
||||||
return 4*np.pi/kw["wavelength"]*np.sin(theta)
|
|
||||||
|
|
||||||
def degToQ(theta,**kw):
|
|
||||||
theta = theta/180.*np.pi
|
|
||||||
return radToQ(theta,**kw)
|
|
||||||
degToQ.__doc__ = radToQ.__doc__
|
|
||||||
|
|
||||||
def qToTheta(q,asDeg=False,**kw):
|
|
||||||
""" Return scattering angle from q (given E or wavelength) """
|
|
||||||
# Energy or wavelength should be in kw
|
|
||||||
assert "E" in kw or "wavelength" in kw
|
|
||||||
# but not both
|
|
||||||
assert not ("E" in kw and "wavelength" in kw)
|
|
||||||
if "E" in kw: kw["wavelength"] = 12.398/kw["E"]
|
|
||||||
theta = np.arcsin(q*kw["wavelength"]/4/np.pi)
|
|
||||||
if asDeg: theta = np.rad2deg(theta)
|
|
||||||
return theta
|
|
||||||
|
|
||||||
def attenuation_length(compound, density=None, natural_density=None,energy=None, wavelength=None):
|
|
||||||
""" extend periodictable.xsf capabilities """
|
|
||||||
import periodictable.xsf
|
|
||||||
if energy is not None: wavelength = periodictable.xsf.xray_wavelength(energy)
|
|
||||||
assert wavelength is not None, "scattering calculation needs energy or wavelength"
|
|
||||||
if (np.isscalar(wavelength)): wavelength=np.array( [wavelength] )
|
|
||||||
n = periodictable.xsf.index_of_refraction(compound=compound,
|
|
||||||
density=density, natural_density=natural_density,
|
|
||||||
wavelength=wavelength)
|
|
||||||
attenuation_length = (wavelength*1e-10)/ (4*np.pi*np.imag(n))
|
|
||||||
return np.abs(attenuation_length)
|
|
||||||
|
|
||||||
def transmission(material='Si',thickness=100e-6, density=None, natural_density=None,energy=None, wavelength=None):
|
|
||||||
""" extend periodictable.xsf capabilities """
|
|
||||||
att_len = attenuation_length(material,density=density,
|
|
||||||
natural_density=natural_density,energy=energy,wavelength=wavelength)
|
|
||||||
return np.exp(-thickness/att_len)
|
|
||||||
|
|
||||||
|
|
||||||
def chargeToPhoton(chargeOrCurrent,material="Si",thickness=100e-6,energy=10,e_hole_pair=3.63):
|
|
||||||
"""
|
|
||||||
Function to convert charge (or current to number of photons (or number
|
|
||||||
of photons per second)
|
|
||||||
|
|
||||||
Parameters
|
|
||||||
----------
|
|
||||||
|
|
||||||
chargeOrCurrent: float or array
|
|
||||||
material : str
|
|
||||||
Used to calculate
|
|
||||||
|
|
||||||
"""
|
|
||||||
# calculate absortption
|
|
||||||
A = 1-transmission(material=material,energy=energy,thickness=thickness)
|
|
||||||
chargeOrCurrent = chargeOrCurrent/A
|
|
||||||
|
|
||||||
e_hole_pair_energy_keV = e_hole_pair*1e-3
|
|
||||||
n_charge_per_photon = energy/e_hole_pair_energy_keV
|
|
||||||
# convert to Q
|
|
||||||
charge_per_photon = n_charge_per_photon*1.60217662e-19
|
|
||||||
nphoton = chargeOrCurrent/charge_per_photon
|
|
||||||
|
|
||||||
if len(nphoton) == 1: nphoton = float(nphoton)
|
|
||||||
return nphoton
|
|
||||||
|
|
||||||
|
|
||||||
def logToScreen():
|
|
||||||
""" 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 (if needed)
|
|
||||||
if len(logging.getLogger('').handlers)==1:
|
|
||||||
logging.getLogger('').addHandler(console)
|
|
||||||
49
xray/xray.py
49
xray/xray.py
|
|
@ -1,49 +0,0 @@
|
||||||
from __future__ import print_function,division
|
|
||||||
import numpy as np
|
|
||||||
import os
|
|
||||||
|
|
||||||
def xrayAttLenght(*args,**kw):
|
|
||||||
from periodictable import xsf
|
|
||||||
n = xsf.index_of_refraction(*args,**kw)
|
|
||||||
if not "wavelength" in kw:
|
|
||||||
wavelength = 12.398/np.asarray(kw["energy"])
|
|
||||||
else:
|
|
||||||
wavelength = np.asarray(kw["wavelength"])
|
|
||||||
attenuation_length = np.abs( (wavelength*1e-10)/ (4*np.pi*np.imag(n)) )
|
|
||||||
return attenuation_length
|
|
||||||
|
|
||||||
|
|
||||||
def xrayFluo(atom,density,energy=7.,length=30.,I0=1e10,\
|
|
||||||
det_radius=1.,det_dist=10.,det_material="Si",det_thick=300,verbose=True):
|
|
||||||
""" compound: anything periodictable would understand
|
|
||||||
density: in mM
|
|
||||||
length: sample length in um
|
|
||||||
energy: in keV, could be array
|
|
||||||
"""
|
|
||||||
import periodictable
|
|
||||||
from periodictable import xsf
|
|
||||||
wavelength = 12.398/energy
|
|
||||||
atom = periodictable.__dict__[atom]
|
|
||||||
# 1e3 -> from mM to M
|
|
||||||
# 1e3 -> from L to cm3
|
|
||||||
# so 1 mM = 1e-3M/L = 1e-6 M/cm3
|
|
||||||
density_g_cm3 = density*1e-6*atom.mass
|
|
||||||
n = xsf.index_of_refraction(atom,density=density_g_cm3,wavelength=wavelength)
|
|
||||||
attenuation_length = xrayAttLenght( atom,density=density_g_cm3,energy=energy )
|
|
||||||
# um-> m: 1e-6
|
|
||||||
fraction_absorbed = 1.-np.exp(-length*1e-6/attenuation_length)
|
|
||||||
if verbose:
|
|
||||||
print("Fraction of x-ray photons absorbed by the sample:",fraction_absorbed)
|
|
||||||
## detector ##
|
|
||||||
det_area = np.pi*det_radius**2
|
|
||||||
det_solid_angle = det_area/(4*np.pi*det_dist**2)
|
|
||||||
if verbose:
|
|
||||||
print("Detector fraction of solid angle:",det_solid_angle)
|
|
||||||
det_att_len = xrayAttLenght(det_material,wavelength=atom.K_alpha)
|
|
||||||
det_abs = 1-np.exp(-det_thick*1e-6/det_att_len)
|
|
||||||
if verbose:
|
|
||||||
print("Detector efficiency (assuming E fluo = E_K_alpha):",det_abs)
|
|
||||||
eff = fraction_absorbed*det_solid_angle*det_abs
|
|
||||||
if verbose:
|
|
||||||
print("Overall intensity (as ratio to incoming one):",eff)
|
|
||||||
return eff
|
|
||||||
Loading…
Reference in New Issue