361 lines
13 KiB
Python
361 lines
13 KiB
Python
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)
|
|
id9.default_extension = '.npz'
|
|
#id9.default_extension = '.h5'
|
|
|
|
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):
|
|
from pyFAI import calibrant
|
|
c120 = calibrant.Cell(a=10.297,b=10.716,c=13.955,alpha=72.22,beta=70.45,gamma=74.13)
|
|
p120 = np.pi*2/ np.asarray( list(c120.d_spacing(dmin=dmin).keys()) ).astype(float)
|
|
p120.sort()
|
|
c300 = calibrant.Cell(a=10.485,b=11.015,c=14.280,alpha=74.61,beta=69.29,gamma=75.29)
|
|
p300 = np.pi*2/ np.asarray( list(c300.d_spacing(dmin=dmin).keys()) ).astype(float)
|
|
p300.sort()
|
|
return xray.storage.DataStorage( dict(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.exists(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"+id9.default_extension
|
|
if not os.path.exists(fname) or force:
|
|
data = azav(folder,doRemoveBaseline=False)
|
|
# 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,
|
|
doRemoveBaseline=True):
|
|
if isinstance(folder,int): folder = g_default_folder%folder
|
|
if isinstance(mask,int):
|
|
files = xray.utils.getFiles(folder,"*.edf*")
|
|
img = xray.azav.pyFAIread(files[0])
|
|
temp = np.ones_like(img,dtype=bool)
|
|
temp[:mask] = False
|
|
mask = temp
|
|
azav = id9.doFolder_azav(folder,nQ=nQ,force=force,mask=mask,saveChi=saveChi)
|
|
if doRemoveBaseline: removeBaseline(folder)
|
|
return azav
|
|
|
|
def datared(folder,monitor=(0.5,4),showPlot=True,errMask=5,chi2Mask=2,
|
|
qlims=(0.5,3),withBaseline=False,forceDatared=False,
|
|
select=slice(None),**kw):
|
|
if isinstance(folder,int): folder = g_default_folder%folder
|
|
def_ext = id9.default_extension
|
|
extension = def_ext if withBaseline else "_nobaseline" + def_ext
|
|
# storage files to use
|
|
diffFile = folder + "/diffs" + extension
|
|
azavFile = folder + "/pyfai_1d" + extension
|
|
|
|
if not os.path.exists(diffFile) or forceDatared:
|
|
data,diffs = id9.doFolder_dataRed(azavFile,monitor=monitor,
|
|
errMask=errMask,chi2Mask=chi2Mask,qlims=qlims,
|
|
outStorageFile=diffFile)
|
|
else:
|
|
data = xray.storage.read(azavFile)
|
|
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):
|
|
if y.ndim == 1: y = (y,)
|
|
ret = dict()
|
|
for yy in y:
|
|
res = xray.peaks.fitPeak(x,yy)
|
|
#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] )
|
|
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")
|
|
data = removeBaseline("../fesalen/fesalen1/run15")
|
|
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:
|
|
log.warn("scaleWithE requires useRatio, setting scaleWithE to False")
|
|
scaleWithE = Fale
|
|
for run in runs:
|
|
t,pos = anaPeakPos(run,plot=False,withBaseline=withBaseline)
|
|
for k,a in zip(pos.keys(),ax):
|
|
label = "run %d %.0f uJ"%(run,runsE[run])
|
|
if useRatio:
|
|
ratio = pos[k].center.value[1:]/pos[k].center.value[0]
|
|
if scaleWithE: ratio = (ratio-1)/runsE[run]+1
|
|
line = a.plot(t[1:],ratio,'-o',label=label)[0]
|
|
else:
|
|
line = a.plot(t[1:],pos[k].center.value[1:],'-o',label=label)[0]
|
|
a.axhline(pos[k].center.value[0],color=line.get_color())
|
|
for interval,a in zip(pos.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}$)")
|