update salen example

xray/example_main_salen.py

@@ -1,6 +1,6 @@
 from __future__ import print_function,division
 import os
-
+import collections
 
 # prepare logging
 import logging
@@ -17,15 +17,28 @@ import mcutils.xray as xray
 from mcutils.xray import id9
 #id9 = xray.id9
 
-from dualtree import dualtree
-
 # use npz files (they can handle more stuff (list of arrays,unicode) than h5py)
 id9.default_extension = '.npz'
 #id9.default_extension = '.h5'
 
 g_default_mask = '../masks/fesalen1_run8_careful.edf'
+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
@@ -35,56 +48,110 @@ def prepareLog():
   formatter = logging.Formatter('%(message)s')
   # tell the handler to use this format
   console.setFormatter(formatter)
-  # add the handler to the root logger
+  # add the handler to the root logger (if needed)
-  logging.getLogger('').addHandler(console)
+  if len(logging.getLogger('').handlers)==1:
+    logging.getLogger('').addHandler(console)
 
 prepareLog()
 
-def findCenter():
+def defineCalibrants(dmin=3):
-  files = xray.utils.getEdfFiles("../fesalen/fesalen1/run8/",nFiles=100)
+  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.
-  xray.azav.find_center(img)
+  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 azav(folder,nQ=1500,force=False,saveChi=True,mask=g_default_mask):
+
+
+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
-  return id9.doFolder_azav(folder,nQ=nQ,force=force,mask=mask,saveChi=saveChi)
+  azav = id9.doFolder_azav(folder,nQ=nQ,force=force,mask=mask,saveChi=saveChi)
-
+  if doRemoveBaseline: removeBaseline(folder)
-
+  return azav
-def removeBaseline(folder,qlims=(0.6,3),max_iter=30):
-  data = azav(folder)
-  idx = (data.q>qlims[0]) & (data.q<qlims[1])
-  data['q'] = data.q[idx]
-  data['data'] = data.data[:,idx]
-  for i in range(len(data.data)):
-    data['data'][i]=data.data[i]-dualtree.baseline(data.data[i],max_iter=max_iter)
-  fname = os.path.splitext(data.filename)[0] + "_nobaseline" + id9.default_extension
-  data.save(fname)
-  return data
 
 def datared(folder,monitor=(0.5,4),showPlot=True,errMask=5,chi2Mask=2,
-            qlims=(0.6,3),withBaseline=False,storageFile='auto',force=False,**kw):
+            qlims=(0.5,3),withBaseline=False,forceDatared=False,
-  # ice contributes to a lot of noise, filter to q<2
+            select=slice(None),**kw):
-  if storageFile == 'auto':
+  if isinstance(folder,int): folder = g_default_folder%folder
-    if withBaseline:
+  def_ext = id9.default_extension
-      storageFile = folder + "/" + "pyfai_1d" + id9.default_extension
+  extension = def_ext if withBaseline else "_nobaseline" + def_ext
-    else:
+  # storage files to use
-      storageFile = folder + "/" + "pyfai_1d_nobaseline" + id9.default_extension
+  diffFile = folder + "/diffs" + extension
-      if not os.path.exists(storageFile) or force: removeBaseline(folder)
+  azavFile = folder + "/pyfai_1d" + extension
-  data = xray.storage.DataStorage(storageFile)
+
-  data,diffs = id9.doFolder_dataRed(folder,storageFile=data,monitor=monitor,
+  if not os.path.exists(diffFile) or forceDatared:
-               errMask=errMask,chi2Mask=chi2Mask,qlims=qlims,**kw)
+    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,
+    xray.utils.plotdiffs(diffs.q,diffs.data,t=diffs.scan,select=select,
-         absSignal=diffs.dataAbsAvAll,absSignalScale=30)
+         absSignal=diffs.dataAbsAvNeg,absSignalScale=10)
     plt.title(folder + " norm %s" % str(monitor))
     plt.figure()
-    xray.utils.plotdiffs(diffs.q,diffs.dataAbsAvScanPoint,t=diffs.scan)
+    xray.utils.plotdiffs(diffs.q,diffs.dataAsAbs,t=diffs.scan,
+         select=select)
     plt.title(folder + " norm %s" % str(monitor))
   return data,diffs
 
@@ -93,14 +160,201 @@ def doall(folder,force=False,removeBaseline=True):
   azav(folder,force=force)
   return datared(folder)
 
-def anaAmplitue(run=6):
+def _getAmplitudeFromSVD(i,comp=0):
-  fname = "../tiox/tiox1/run%d/diffs.npz" % run
+  u,s,v = np.linalg.svd(i)
-  data  = xray.storage.DataStorage(fname)
+  return u[:,comp]*s[comp]
-  ranges = ( (1.75,1.85), (2.2,2.4), (3.25,3.4) )
+
+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)
-  fig,ax = plt.subplots(nPlot,1,sharex=True)
+  if plot: fig,ax = plt.subplots(nPlot,1,sharex=True)
-  for r,a in zip(ranges,ax):
+  amp = collections.OrderedDict()
-    idx = (data.q>r[0]) & (data.q<r[1])
+  for num,r in enumerate(ranges):
-    amplitude = np.abs(data.data[:,idx]).mean(axis=1)
+    idx = (diffs.q>r[0]) & (diffs.q<r[1])
-    a.plot(data.scan,amplitude,'-o')
+    #amplitude = _getAmplitudeFromLinFit( diffs.data[:,idx] )
-    a.set_title("Range %s"%(str(r)))
+    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}$)")