first git.ipr commit

.gitignore

@@ -1,3 +1,6 @@
+mecl3616_output/
+xppl3716_output/
+
 # ---> Python
 # Byte-compiled / optimized / DLL files
 __pycache__/

.ipynb_checkpoints/HOW-TO-checkpoint.ipynb

@@ -0,0 +1,6 @@
+{
+ "cells": [],
+ "metadata": {},
+ "nbformat": 4,
+ "nbformat_minor": 0
+}

NOTES.txt

@@ -0,0 +1,6 @@
+run 13, fit on first shot, foms:
+{'fom': array([ 0.05494294,  0.06082832,  0.0530633 ,  0.04097997,  0.04311072,
+         0.04868919,  0.03658052,  0.05205008,  0.02965361,  0.0349829 ,
+         0.05938427,  0.04889882,  0.04420127,  0.05612349,  0.05459789,
+         0.04026893,  0.06067866,  0.04404101,  0.05147894,  0.06231822]),
+

alignment.py

@@ -0,0 +1,589 @@
+from __future__ import print_function,division
+from skimage import transform as tf
+from scipy.ndimage import gaussian_filter1d as g1d
+import mcutils as mc
+import joblib
+import collections
+import numpy as np
+import matplotlib.pyplot as plt
+import time
+
+# /--------\
+# |        |
+# | UTILS  |
+# |        |
+# \--------/
+
+#defaultE = np.arange(1024)*0.12+7060
+defaultE = (np.arange(1024)-512)*0.189+7123
+
+__i = np.arange(1024)
+__x = (__i-512)/512
+
+g_fit_default_kw = dict(
+    intensity        = 1.,
+    error_intensity  = 0.02,
+    transx           = 0,
+    error_transx     = 3,
+    limit_transx     = ( -400,400 ),
+    transy           = 0,
+    error_transy     = 3,
+    limit_transy     = ( -50,50 ),
+    rotation         = 0.01,
+    error_rotation   = 0.005,
+    limit_rotation   = (-0.06,0.06),
+    scalex           = 1,
+    limit_scalex     = (0.4,1.2),
+    error_scalex     = 0.05,
+    scaley           = 1,
+    error_scaley     = 0.05,
+    limit_scaley     = (0.8,1.2),
+    shear            = 0.01,
+    error_shear      = 0.001,
+    limit_shear      = (-0.2,0.2),
+    igauss1cen       = 512,
+    error_igauss1cen = 2.,
+    fix_igauss1cen   = True,
+    igauss1sig       = 4000.,
+    error_igauss1sig = 2.,
+    fix_igauss1sig   = True,
+    igauss2cen       = 512,
+    error_igauss2cen = 2.,
+    fix_igauss2cen   = True,
+    igauss2sig       = 4000.,
+    error_igauss2sig = 2.,
+    fix_igauss2sig   = True,
+    iblur1           = 0,
+    limit_iblur1     = (0,20),
+    error_iblur1     = 0.02,
+    fix_iblur1       = True
+)
+
+
+def rebin1D(a,shape):
+  n0 = a.shape[0]//shape
+  sh = shape,n0
+  return a[:n0*shape].reshape(sh).mean(1)
+
+
+kw_2dplot = dict(
+  interpolation = "none",
+  aspect        = "auto",
+  cmap          = plt.cm.viridis
+  )
+
+fit_ret = collections.namedtuple("fit_ret",["fit_result","init_pars","final_pars",\
+  "final_transform1","final_transform2","im1","im2","p1","p2","fom","ratio","tneeded"] )
+
+def calcFOM(p1,p2,ratio):
+  idx = ( p1>p1.max()/10 ) & (p2>p2.max()/10)
+  ratio = ratio[idx]
+  return ratio.std()/ratio.mean()
+  
+def subtractBkg(imgs,nPix=100,bkg_type="line"):
+  if imgs.ndim == 2: imgs = imgs[np.newaxis,:]
+  imgs = imgs.astype(np.float)
+  if bkg_type == "line":
+    bkg = np.median(imgs[:,:nPix,:],axis=1)
+    imgs = imgs-bkg[:,np.newaxis,:]
+  elif bkg_type == "corner":
+    q1 = imgs[:,:nPix,:nPix].mean(-1).mean(-1)
+    imgs[:,:512,:512]-=q1[:,np.newaxis,np.newaxis]
+    q2 = imgs[:,:nPix,-nPix:].mean(-1).mean(-1)
+    imgs[:,:512,-512:]-=q2[:,np.newaxis,np.newaxis]
+    q3 = imgs[:,-nPix:,-nPix:].mean(-1).mean(-1)
+    imgs[:,-512:,-512:]-=q3[:,np.newaxis,np.newaxis]
+    q4 = imgs[:,-nPix:,:nPix].mean(-1).mean(-1)
+    imgs[:,-512:,:512]-=q4[:,np.newaxis,np.newaxis]
+  elif bkg_type is None:
+    if imgs.ndim == 2: imgs = imgs[np.newaxis,:].astype(np.float)
+  else:
+    print("Background subtraction '%s' Not impleted"%bkg_type)
+  return imgs
+
+def getCenterOfMax(img,axis=0,threshold=0.05):
+  img = img.copy()
+  img[img<img.max()*threshold] = 0
+  if axis == 1: img=img.T
+  p   = img.mean(1)
+  x   = np.arange(img.shape[0])
+  return np.sum(x*p)/np.sum(p)
+
+
+def findRoi(img,height=100,axis=0):
+  c = int( getCenterOfMax(img,axis=axis) )
+  roi = slice(c-height//2,c+height//2)
+  return roi
+
+# /--------------------\
+# |                    |
+# |     PLOTS & CO.    |
+# |                    |
+# \--------------------/
+
+def plotShot(im1,im2,transf1=None,transf2=None,fig=None,ax=None,res=None,E=defaultE,save=None):
+  if transf1 is not None: im1 = transf1.transformImage(im1)
+  if transf2 is not None: im2 = transf2.transformImage(im2)
+  if fig is None and ax is None:
+    fig = plt.subplots(2,3,figsize=[7,5],sharex=True)[0]
+    ax = fig.axes
+  elif fig is not None:
+    ax = fig.axes
+  if E is None: E=np.arange(im1.shape[1])
+  n = im1.shape[0]
+  ax[0].imshow(im1,**kw_2dplot,extent=(E[0],E[-1],0,n))
+  ax[1].imshow(im2,**kw_2dplot,extent=(E[0],E[-1],0,n))
+  ax[2].imshow(im1-im2,**kw_2dplot,extent=(E[0],E[-1],0,n))
+  if res is None:
+    p1 = np.nansum(im1,axis=0)
+    p2 = np.nansum(im2,axis=0)
+    pr = p1/p2
+  else:
+    p1 = res.p1; p2 = res.p2; pr = res.ratio
+  ax[3].plot(E,p1,lw=3)
+  ax[4].plot(E,p1,lw=1)
+  ax[4].plot(E,p2,lw=3)
+  idx = (p1>p1.max()/10.) & (p2>p2.max()/10)
+  ax[5].plot(E[idx],pr[idx])
+  if res is not None:
+    ax[5].set_title("FOM: %.2f"%res.fom)
+  else:
+    ax[5].set_title("FOM: %.2f"% calcFOM(p1,p2,pr))
+  if (save is not None) and (save is not False): plt.savefig(save,transparent=True,dpi=500)
+  return fig
+
+def plotRatios(r,fig=None,E=defaultE,save=None):
+  if fig is None: fig = plt.subplots(2,1,sharex=True)[0]
+  ax = fig.axes
+  n = r.shape[0]
+  i = ax[0].imshow(r,**kw_2dplot,extent=(E[0],E[-1],0,n))
+  i.set_clim(0,1.2)
+  idx = np.random.random_integers(0,n-1)
+  ax[1].plot(E,r[idx],label="Shot n %d"%idx)
+  ax[1].plot(E,np.nanmedian(r[:10],axis=0),label="median 10 shots")
+  ax[1].plot(E,np.nanmedian(r,axis=0),label="all shots")
+  ax[1].legend()
+  ax[1].set_ylim(0,1.5)
+  ax[1].set_xlabel("Energy")
+  ax[1].set_ylabel("Transmission")
+  ax[0].set_ylabel("Shot num")
+  if (save is not None) and (save is not False): plt.savefig(save,transparent=True,dpi=500)
+
+def plotSingleShots(r,nShots=10,fig=None,E=defaultE,save=None,ErangeForStd=(7090,7150)):
+  if fig is None: fig = plt.subplots(2,1,sharex=True)[0]
+  ax = fig.axes
+  for i in range(nShots):
+    ax[0].plot(E,r[i]+i)
+  ax[0].set_ylim(0,nShots+0.5)
+  av = (1,3,10,30,100)
+  good = np.nanmedian(r,0)
+  for i,a in enumerate(av):
+   m = np.nanmedian(r[:a],0)
+   idx = (E>ErangeForStd[0]) & (E<ErangeForStd[1])
+   fom = np.nanstd( m[idx]/good[idx] ) 
+   print("n shots %d, std %.2f"%(a,fom) )
+   ax[1].plot(E,m+i,label="%d shots, std :%.2f"%(a,fom))
+  ax[1].legend()
+  ax[1].set_ylim(0,len(av)+0.5)
+  ax[1].set_xlabel("Energy")
+  ax[1].set_ylabel("Transmission")
+  if (save is not None) and (save is not False): plt.savefig(save,transparent=True,dpi=500)
+
+
+# /--------------------\
+# |                    |
+# |  TRANSFORM & CO.   |
+# |                    |
+# \--------------------/
+
+
+def transformImage(img,transform=None,iblur=None,intensity=1,\
+    igauss=None,orderRotation=1,show=False):
+  """ Transform is the geometrical (affine) transform)
+      blur is a tuple (blurx,blury); if a single value used for energy axis
+      i    is to correct the itensity
+      igauss is to multiply the intensity with a gaussian along the energy axis (=0); if a single assumed centered at center of image, or has to be a tuple (cen,witdth)
+  """
+  if transform is not None and not np.all(transform.params==np.eye(3)) :
+    try:
+      t = np.linalg.inv(transform.params)
+      i = tf._warps_cy._warp_fast(img,t,order=orderRotation)
+    except np.linalg.LinAlgError:
+      print("Image transformation failed, returning original image")
+      i = img.copy()
+  else:
+    i = img.copy()
+  i *= intensity
+  if iblur is not None:
+    if isinstance(iblur,(int,float)):
+      iblur = (iblur,None)
+    if (iblur[0] is not None) and (iblur[0]>0): i = g1d(i,iblur[0],axis=1)
+    if (iblur[1] is not None) and (iblur[1]>0): i = g1d(i,iblur[1],axis=0)
+  if igauss is not None:
+    if isinstance(igauss,(int,float)):
+      igauss = (__i[-1]/2,igauss); # if one parameter only, assume it is centered on image
+    g = mc.gaussian(__i,x0=igauss[0],sig=igauss[1],normalize=False)
+    i *= g
+  if show: plotShot(img,i)
+  return i
+
+
+class SpecrometerTransformation(object):
+  def __init__(self,translation=(0,0),scale=(1,1),rotation=0,shear=0,
+               intensity=1,igauss=None,iblur=None):
+    self.affineTransform = getTransform(translation=translation,
+      scale = scale,rotation=rotation,shear=shear)
+    self.intensity=intensity
+    self.igauss = igauss
+    self.iblur   = iblur
+
+  def update(self,**kw):
+    # current transformation, necessary because skimage transformation do not support
+    # setting of attributes
+    names = ["translation","scale","rotation","shear"]
+    trans_dict = dict( [(n,getattr(self.affineTransform,n)) for n in names] )
+    for k,v in kw.items():
+      if hasattr(self,k):
+        setattr(self,k,v)
+      elif k in names:
+        trans_dict[k] = v
+    self.affineTransform = getTransform(**trans_dict)
+
+  def transformImage(self,img,orderRotation=1,show=False):
+    return transformImage(img,self.affineTransform,iblur=self.iblur,\
+      intensity=self.intensity,igauss=self.igauss,orderRotation=orderRotation,\
+      show=show)
+
+def saveAlignment(fname,transform,roi1,roi2):
+  np.save(fname, dict( transform = transform, roi1=roi1, roi2 = roi2 ) )
+
+def loadAlignment(fname):
+  return np.load(fname).item()
+
+def unravel_results(res,getBest=False):
+  out  = dict()
+  parnames = res[0].fit_result.parameters
+  out["parameters"] = dict()
+  for n in parnames:
+    out["parameters"][n] = np.asarray( [r.final_pars[n] for r in res])
+  out["ratio"] = np.asarray( [r.ratio for r in res])
+  out["p1"] =    np.asarray( [r.p1 for r in res]  )
+  out["p2"] =    np.asarray( [r.p2 for r in res]  )
+  out["fom"] =   np.asarray( [r.fom for r in res] )
+  return out
+
+def getTransform(translation=(0,0),scale=(1,1),rotation=0,shear=0):
+  t= tf.AffineTransform(scale=scale,rotation=rotation,shear=shear,\
+         translation=translation)
+  return t
+
+def findTransform(p1,p2,ttype="affine"):
+  return tf.estimate_transform(ttype,p1,p2)
+
+
+#__i = np.arange(2**12)/2**11
+#def _transformToIntensitylDependence(transform):
+#  c = 1.
+#  if hasattr(transform,"i_a"):
+#    c += transform.i_a*_i
+#  return c
+
+
+
+# /--------\
+# |        |
+# |   FIT  |
+# |        |
+# \--------/
+
+
+def transformIminuit(im1,im2,init_transform=dict(),show=False,verbose=True,zeroThreshold=0.05,doFit=True):
+  import iminuit
+  assert im1.dtype == im2.dtype
+  t0 = time.time()
+  # create local copy we can mess up with
+  im1_toFit = im1.copy() 
+  im2_toFit = im2.copy() 
+
+  # set anything below the 5% of the max to zero (one of the two opal is noisy)
+  im1_toFit[im1_toFit<im1.max()*zeroThreshold] = 0
+  im2_toFit[im2_toFit<im2.max()*zeroThreshold] = 0
+  p1 = im1.mean(0)
+  p2 = im2.mean(0)
+
+  # set errorbar
+  err = 3.
+
+  def transforms(intensity,
+    igauss1cen,igauss1sig,iblur1,
+    scalex,scaley,rotation,transx,transy,shear,
+    igauss2cen,igauss2sig):
+    t1 = SpecrometerTransformation(translation=(transx,transy),scale=(scalex,scaley),\
+         rotation=rotation,shear=shear,intensity=intensity,igauss=(igauss1cen,igauss1sig),\
+         iblur=iblur1)
+    t2 = SpecrometerTransformation(igauss=(igauss2cen,igauss2sig))
+    return t1,t2
+
+  def model(intensity,
+    igauss1cen,igauss1sig,iblur1,
+    scalex,scaley,rotation,transx,transy,shear,
+    igauss2cen,igauss2sig):
+
+    t1,t2 = transforms(intensity,
+    igauss1cen,igauss1sig,iblur1,
+    scalex,scaley,rotation,transx,transy,shear,
+    igauss2cen,igauss2sig)
+    return t1.transformImage(im1_toFit),t2.transformImage(im2_toFit)
+  
+  def chi2(intensity,
+    igauss1cen,igauss1sig,iblur1,
+    scalex,scaley,rotation,transx,transy,shear,
+    igauss2cen,igauss2sig):
+
+    i1,i2 = model(intensity,     \
+    igauss1cen,igauss1sig,iblur1, \
+    scalex,scaley,rotation,transx,transy,shear, \
+    igauss2cen,igauss2sig)
+    d = (i1-i2)/err
+    return np.sum(d*d)
+
+  # set default initial stepsize and limits
+  r = im2.mean(0).sum()/im1.mean(0).sum()
+
+  default_kw = g_fit_default_kw.copy()
+  default_kw["intensity"] = r
+
+  init_kw = dict()
+  if isinstance(init_transform,dict):
+    init_kw = init_transform.copy()
+  elif isinstance(init_transform,iminuit._libiminuit.Minuit):
+    init_kw = init_transform.fitarg.copy()
+  elif isinstance(init_transform,tf.AffineTransform):
+    init_kw["transx"],init_kw["transy"] = init_transform.translation
+    init_kw["scalex"],init_kw["scaley"] = init_transform.scale
+    init_kw["rotation"] = init_transform.rotation
+    init_kw["shear"] = init_transform.shear
+  if "intensity" in init_kw and init_kw["intensity"] == "auto":
+    r = im2.mean(0).sum()/im1.mean(0).sum()
+    init_kw["intensity"]= r
+
+  kw = default_kw.copy()
+  kw.update(init_kw)
+  kw["fix_shear"]=True
+
+
+  tofix = ("scalex","scaley","rotation","shear")
+  kw_tofix = dict( [("fix_%s"%p,True) for p in tofix] )
+  kw.update(kw_tofix)
+  imin  = iminuit.Minuit(chi2,errordef=1.,**kw)
+  imin.set_strategy(1)
+  init_params = imin.fitarg.copy()
+
+  if show:
+    i1,i2 = model(*imin.args)
+    plotShot(i1,i2)
+    fig = plt.gcf()
+    fig.text(0.5,0.9,"Initial Pars")
+    input("Enter to start fit")
+
+  if doFit: imin.migrad()
+
+  pars = imin.fitarg.copy()
+  kw_tofree = dict( [("fix_%s"%p,False) for p in tofix] )
+  pars.update(kw_tofree)
+  imin  = iminuit.Minuit(chi2,errordef=1.,**pars)
+  imin.set_strategy(1)
+
+  if doFit: imin.migrad()
+
+  final_params = imin.fitarg.copy()
+  t1,t2 = transforms(*imin.args)
+  i1,i2 = model(*imin.args)
+  if show:
+    plotShot(i1,i2)
+    fig = plt.gcf()
+    fig.text(0.5,0.9,"Final Pars")
+
+  p1 = np.nansum(i1,axis=0)
+  p2 = np.nansum(i2,axis=0)
+  r  = p2/p1
+  idx = p1>np.nanmax(p1)/10.
+  fom = calcFOM(p1,p2,r)
+  return fit_ret(
+    fit_result   = imin,
+    init_pars    = init_params,
+    final_pars    = final_params,
+    final_transform1 = t1,
+    final_transform2 = t2,
+    im1          = i1,
+    im2          = i2,
+    p1           = p1,
+    p2           = p2,
+    ratio        = r,
+    fom          = fom,
+    tneeded      = time.time()-t0
+  )
+
+ 
+# /--------\
+# |        |
+# | MANUAL |
+# | ALIGN  |
+# |        |
+# \--------/
+class GuiAlignment(object):
+  def __init__(self,im1,im2,autostart=False):
+    self.im1 = im1
+    self.im2 = im2
+    self.f,self.ax=plt.subplots(1,2)
+    self.ax[0].imshow(im1,aspect="auto")
+    self.ax[1].imshow(im2,aspect="auto")
+    self.transform = None
+    if autostart:
+      return self.start()
+    else:
+      print("Zoom first then use the .start method")
+
+  def OnClick(self,event):
+    if event.button == 1:
+      if self._count % 2 == 0:
+        self.im1_p.append( (event.xdata,event.ydata) )
+        print("Added",event.xdata,event.ydata,"to im1")
+      else:
+        self.im2_p.append( (event.xdata,event.ydata) )
+        print("Added",event.xdata,event.ydata,"to im2")
+      self._count += 1
+    elif event.button == 2:
+      # neglect middle click
+      return
+    elif event.button == 3:
+      self.done = True
+      return
+    else:
+      return
+
+  def start(self):
+    self._nP = 0
+    self._count = 0
+    self.im1_p = []
+    self.im2_p = []
+    self.done = False
+    cid_up = self.f.canvas.mpl_connect('button_press_event', self.OnClick)
+    print("Press right button to finish")
+    while not self.done:
+      if self._count % 2 == 0:
+        print("Select point %d for left image"%self._nP)
+      else:
+        print("Select point %d for right image"%self._nP)
+      self._nP = self._count //2
+      plt.waitforbuttonpress()
+    self.im1_p = np.asarray(self.im1_p)
+    self.im2_p = np.asarray(self.im2_p)
+    self.transform = findTransform(self.im1_p,self.im2_p)
+    self.transform.intensity = 1.
+    return self.transform
+
+  def show(self):
+    if self.transform is None:
+      print("Do the alignment first (with .start()")
+      return
+    # just to save some tipying
+    im1 = self.im1; im2 = self.im2
+    im1_new = transformImage(im1,self.transform)
+    f,ax=plt.subplots(1,3)
+    ax[0].imshow(im1,aspect="auto")
+    ax[1].imshow(im1_new,aspect="auto")
+    ax[2].imshow(im2,aspect="auto")
+
+  def save(self,fname):
+    if self.transform is None:
+      print("Do the alignment first (with .start()")
+      return
+    else:
+      np.save(fname,self.transform)
+
+  def load(self,fname):
+    self.transform = np.load(fname).item()
+
+
+def getAverageTransformation(out):
+  res = unravel_results(out)
+  # get average parameters
+  tx = np.median(res["transx"])
+  ty = np.median(res["transy"])
+  sx = np.median(res["scalex"])
+  sy = np.median(res["scaley"])
+  r  = np.median(res["rotation"])
+  sh = np.median(res["shear"])
+  inten  = np.median(res["intensity"])
+  t  = getTransform( (tx,ty),(sx,sy),r,sh,intensity=inten )
+  return t
+
+def checkAverageTransformation(out,imgs1):
+  t,inten = getAverageTransformation(out)
+  res["ratio_av"] = []
+  for shot,values in out.items():
+    i = transformImage(imgs1[shot],t)
+    p1 = np.nansum(i,axis=0)
+    r  = p1/values.p2
+    res["ratio_av"].append( r )
+  res["ratio_av"] = np.asarray(res["ratio_av"])
+  return res
+ 
+
+g_lastpars = None
+
+def clearCache():
+  globals()["g_lastpars"]=None
+
+def doShot(i1,i2,init_pars,doFit=True,show=False):
+  #if g_lastpars is not None: init_pars = g_lastpars
+  r = transformIminuit(i1,i2,init_pars,show=show,verbose=False,doFit=doFit)
+  return r
+
+
+def doShots(imgs1,imgs2,initpars,nJobs=16,doFit=False,returnBestTransform=False):
+  clearCache()
+  N = imgs1.shape[0]
+  pool = joblib.Parallel(backend="threading",n_jobs=nJobs) \
+    (joblib.delayed(doShot)(imgs1[i],imgs2[i],initpars,doFit=doFit) for i in range(N))
+  if returnBestTransform:
+    idx = np.argmin( np.abs([p.fom for p in pool]) ); # abs because sometime fit screws up the give negative spectra...
+    print("FOM for best alignment %.2f"%pool[idx].fom)
+    return pool,pool[idx].final_transform1
+  else:
+    return pool
+#  out = collections.OrderedDict( enumerate(pool) )
+#  return out
+    
+# /--------\
+# |        |
+# | TEST   |
+# | ALIGN  |
+# |        |
+# \--------/
+
+
+
+def testDiling(N=100,roi=slice(350,680),doGUIalignment=False,nJobs=4,useIminuit=True):
+  import xppll37_mc
+  globals()["g_lastpars"]=None
+  d = xppll37_mc.readDilingDataset()
+  im1 = xppll37_mc.subtractBkg(d.opal1[:N])[:,roi,:]
+  im2 = xppll37_mc.subtractBkg(d.opal2[:N])[:,roi,:]
+  N   = im1.shape[0]; # redefie N in case it reads less than N
+  # to manual alignment
+  if doGUIalignment:
+    a = GuiAlignment(im1[0],im2[0])
+    input("Ok to start")
+    init_pars = a.start()
+    np.save("gui_align_transform.npy",init_pars)
+  else:
+    init_pars = np.load("gui_align_transform.npy").item()
+  pool = joblib.Parallel(backend="threading",n_jobs=nJobs,verbose=20) \
+    (joblib.delayed(doShot)(im1[i],im2[i],init_pars,useIminuit=useIminuit) for i in range(N))
+  out = collections.OrderedDict( enumerate(pool) )
+  return out
+
+if __name__ == '__main__':
+  pass

make_little_data.py

@@ -0,0 +1,16 @@
+import os
+import numpy as np
+import xppll37_mc
+import sys
+
+
+def do(run,N):
+  r = xppll37_mc.readDataset(run)
+  o1 = r.opal0[:N]
+  o2 = r.opal1[:N]
+  np.savez("littleData/run%04d.npz" % run, opal0 = o1, opal1 = o2 ) 
+
+if __name__ == "__main__":
+  run = int(sys.argv[1])
+  N=100
+  do(run,N)

online/sMhelper.py

@@ -0,0 +1,276 @@
+import sys
+import zmq
+import numpy as np
+import time
+import pickle
+import alignment
+import matplotlib.pyplot as plt
+import threading
+import datetime
+import copy
+
+
+
+
+def histVec(v,oversample=1):
+  v = np.atleast_1d(v)
+  v = np.unique(v)
+  vd = np.diff(v)
+  vd = np.hstack([vd[0],vd])
+  #vv = np.hstack([v-vd/2,v[-1]+vd[-1]/2])
+  vv = np.hstack([v-vd/2.,v[-1]+vd[-1]/2.])
+  if oversample>1:
+    vvo = []
+    for i in range(len(vv)-1):
+      vvo.append(np.linspace(vv[i],vv[i+1],oversample+1)[:-1])
+    vvo.append(vv[-1])
+    vv = np.array(np.hstack(vvo))
+  return vv
+
+def subtractBkg(imgs,nPix=100,dKtype='corners'):
+  """ Opals tend to have different backgroud for every quadrant """
+  if dKtype is 'corners':
+    if imgs.ndim == 2: imgs = imgs[np.newaxis,:]
+    imgs = imgs.astype(np.float)
+    q1 = imgs[:,:nPix,:nPix].mean(-1).mean(-1)
+    imgs[:,:512,:512]-=q1[:,np.newaxis,np.newaxis]
+    q2 = imgs[:,:nPix,-nPix:].mean(-1).mean(-1)
+    imgs[:,:512,-512:]-=q2[:,np.newaxis,np.newaxis]
+    q3 = imgs[:,-nPix:,-nPix:].mean(-1).mean(-1)
+    imgs[:,-512:,-512:]-=q3[:,np.newaxis,np.newaxis]
+    q4 = imgs[:,-nPix:,:nPix].mean(-1).mean(-1)
+    imgs[:,-512:,:512]-=q4[:,np.newaxis,np.newaxis]
+  elif dKtype is 'stripes':
+    if imgs.ndim == 2: imgs = imgs[np.newaxis,:]
+    imgs = imgs.astype(np.float)
+    s1 =  imgs[:,:nPix,:].mean(-2)
+
+
+
+
+  return np.squeeze(imgs)
+
+def getData():
+  t0 = time.time()
+  context = zmq.Context()
+  socket = context.socket(zmq.REQ)
+  socket.connect('tcp://daq-xpp-mon06:12322')
+  #socket.setsockopt(zmq.SUBSCRIBE, b'')
+  while True:
+    socket.send(b"Request")
+    ret    = socket.recv()
+    ret = pickle.loads(ret, encoding='latin1')
+    print('received',ret.keys(),time.time()-t0)
+    t0 = time.time()
+
+if __name__ == "__main__":
+  getData()
+
+class SpecHandler(object):
+  def __init__(self,connectString='tcp://daq-xpp-mon06:12322',spec1name='opal0',spec2name='opal1',roi=[0,1024,0,1024]):
+    self.connectString = connectString
+    self.spec1name = spec1name
+    self.spec2name = spec2name
+    self.surveyplot = Surveyplot(spec1name=spec1name,spec2name=spec2name)
+    self.roi = roi
+    self.projsimple = ProjSimple(spec1name=spec1name,spec2name=spec2name)
+    self._rawContinuousTime = None
+    self.lastDat = None
+    self.openSocket()
+    self.dataCollector = []
+    self.runningPlot = RunningPlot(self.dataCollector)
+
+
+  def openSocket(self):
+    context = zmq.Context()
+    self.context = context
+    self.socket = context.socket(zmq.REQ)
+    self.socket.connect(self.connectString)
+    #self.socket.setsockopt(zmq.SUBSCRIBE, b'')
+
+  def closeSocket(self):
+    del self.context
+    del self.socket
+
+  def getData(self):
+    self.socket.send(b"Request")
+    ret    = self.socket.recv()
+    ret = pickle.loads(ret, encoding='latin1')
+    for sn in [self.spec1name,self.spec2name]:
+      ret[sn] = np.squeeze(alignment.subtractBkg(ret[sn], nPix=100, bkg_type='line'))
+    self.lastDat = ret
+    return ret
+
+  def getRaw(self,repoenConnection=False,doAlign=False,show=False,doFit=False):
+    if doFit is True: doFit='iminuit'
+
+
+    if repoenConnection:
+      self.closeSocket()
+      self.openSocket()
+    dat = self.getData()
+    im1 = dat[self.spec1name]; im2 = dat[self.spec2name]
+    if doAlign:
+      #t = np.load("gui_align_transform_xppl3716.npy").item()
+      if hasattr(self,'transformer'):
+        algn = self.transformer
+      else:
+        algn = alignment.loadAlignment('last_trafo.npy')
+      t = algn['transform']
+      roi1 = algn['roi1']
+      roi2 = algn['roi2']
+      
+      r = alignment.doShot( im1[roi1,:],im2[roi2,:],t, show=show, doFit=doFit)
+      self.transformer = dict(transform=r.final_transform,roi1=roi1,roi2=roi2)
+      alignment.saveAlignment('last_trafo.npy',r.final_transform,roi1,roi2)
+
+      im1 = r.im1; im2 = r.im2
+      showDiff = True
+      showRatio = True
+    else:
+      showDiff = False
+      showRatio = False
+    self.surveyplot.plotImages(im1,im2,showDiff=showDiff,showRatio=showRatio)
+    self.projsimple.plotProfiles(im1,im2)
+    if doAlign:
+      thres = 0.05
+      im1[im1<thres*np.max(im1.ravel())] = np.nan
+      im2[im2<thres*np.max(im2.ravel())] = np.nan
+      self.dataCollector.append(\
+	  dict( time=datetime.datetime.now(),
+	    fom=r.fom, 
+	    ratProj = np.nansum(im2/im1,axis=0),
+	    im1Proj = np.nansum(im1,axis=0),
+	    im2Proj = np.nansum(im2,axis=0)))
+    self.runningPlot.updatePlot()
+      
+
+  def alignFeatures(self):
+    im1 = copy.copy(self.lastDat[self.spec1name])
+    im2 = copy.copy(self.lastDat[self.spec2name])
+    roi1 = alignment.findRoi(im1)
+    roi2 = alignment.findRoi(im2)
+
+    tra = alignment.GuiAlignment(im1[roi1,:],im2[roi2,:],autostart=False)
+    self.transformer = dict(transform=tra.start(),roi1=roi1,roi2=roi2)
+
+  
+  def getRawContinuuous(self,sleepTime):
+    self._rawContinuousTime = sleepTime
+    if not hasattr(self,'_rawContinuousThread'):
+      def upd():
+        while not self._rawContinuousTime is None:
+          self.getRaw()
+          plt.draw()
+          time.sleep(self._rawContinuousTime)
+
+      self._rawContinuousThread = threading.Thread(target=upd)
+      self._rawContinuousThread.start()
+
+
+
+class Surveyplot(object):
+  def __init__(self,spec1name='spec1',spec2name='spec2'):
+    self.fig,self.axs = plt.subplots(4,1,sharex=True,sharey=True)
+    self.axs[0].set_title(spec1name)
+    self.axs[1].set_title(spec2name)
+    self.axs[2].set_title("Difference")
+    self.axs[3].set_title("Ratio")
+
+  def plotImages(self,img1,img2,showDiff=False,showRatio=False):
+    if hasattr(self,'i1'):
+      self.i1.set_data(img1)
+    else:
+      self.i1 = self.axs[0].imshow(img1,origin='lower',interpolation='none')
+    if hasattr(self,'i2'):
+      self.i2.set_data(img2)
+    else:
+      self.i2 = self.axs[1].imshow(img2,origin='lower',interpolation='none')
+    if showDiff:
+      tdiff = img2-img1
+      if hasattr(self,'idiff'):
+        self.idiff.set_data(tdiff)
+      else:
+        self.idiff = self.axs[2].imshow(tdiff,interpolation='none',origin='lower')
+        lms = np.percentile(tdiff,[30,70])
+        self.idiff.set_clim(lms)
+    if showRatio:
+      tratio = img2/img1
+      if hasattr(self,'iratio'):
+        self.iratio.set_data(tratio)
+      else:
+        self.iratio = self.axs[3].imshow(tratio,interpolation='none',origin='lower')
+        lms = np.percentile(tratio,[30,70])
+        self.iratio.set_clim(lms)
+
+class ProjSimple(object):
+  def __init__(self,spec1name='spec1',spec2name='spec2',roi1=[0,1024,0,1024],roi2=[0,1024,0,1024]):
+    self.fig,self.axs = plt.subplots(2,1)
+    self.roi1 = roi1
+    self.roi2 = roi2
+    self.spec1name=spec1name
+    self.spec2name=spec2name
+  
+  #def getROI(self,specNo):
+  def _roiit(self,img,roi):
+    return img[roi[0]:roi[1],roi[2]:roi[3]]
+  
+  def plotProfiles(self,img1,img2):
+    prof1 = np.nansum(self._roiit(img1,self.roi1),0)
+    prof2 = np.nansum(self._roiit(img2,self.roi2),0)
+    if hasattr(self,'l1'):
+      self.l1.set_ydata(prof1)
+    else:
+      self.l1 = self.axs[0].plot(prof1,label=self.spec1name)[0]
+    if hasattr(self,'l2'):
+      self.l2.set_ydata(prof2)
+    else:
+      self.l2 = self.axs[0].plot(prof2,label=self.spec2name)[0]
+      plt.legend()
+
+    if hasattr(self,'lrat'):
+      self.lrat.set_ydata(prof2/prof1)
+    else:
+      self.lrat = self.axs[1].plot(prof2/prof1,label='ratio')[0]
+      self.axs[1].set_ylim(0,2)
+
+
+class RunningPlot(object):
+  def __init__(self,dataCollector):
+    self.dataCollector = dataCollector
+    self.fig,self.axs = plt.subplots(1,1)
+    
+  def updatePlot(self):
+    if len(self.dataCollector)>0:
+      times = np.asarray(([i['time'] for i in self.dataCollector]))
+      foms = np.asarray(([i['fom'] for i in self.dataCollector]))
+      im1Proj = np.asarray(([i['im1Proj'] for i in self.dataCollector]))
+      im2Proj = np.asarray(([i['im2Proj'] for i in self.dataCollector]))
+
+      if hasattr(self,'fomline'):
+        self.fomline.set_xdata(times)
+        self.fomline.set_ydata(foms)
+	#self.axs[0].autoscale(enable=True,axis='x')
+
+      else:
+        self.fomline = self.axs.plot(times,foms,'o-')[0]
+      #if hasattr(self,'ratioimg'):
+	#self.ratioimg.set_data()
+	#self.ratioimg.set_ydata(foms)
+      #else:
+	#self.axs[0].plot(times,foms,'o-')
+
+
+   
+
+  
+  #def plotOrUpdate(img1,img2):
+    #if hasattr(self,i1):
+      #self.i1.set_data(img1)
+    #else:
+      #self.i1 = self.axs.imshow(img1,interpolate='none',origin='bottom')
+    #if hasattr(self,i1):
+      #self.i1.set_data(img1)
+    #else:
+      #self.i1 = self.axs.imshow(img1,interpolate='none',origin='bottom')
+

online/sMhelper_dev.py

@@ -0,0 +1,297 @@
+import sys
+import zmq
+import numpy as np
+import time
+import pickle
+import alignment
+import matplotlib.pyplot as plt
+import threading
+import datetime
+import copy
+
+plt.rcParams['image.cmap'] = 'viridis'
+
+def histVec(v,oversample=1):
+  v = np.atleast_1d(v)
+  v = np.unique(v)
+  vd = np.diff(v)
+  vd = np.hstack([vd[0],vd])
+  #vv = np.hstack([v-vd/2,v[-1]+vd[-1]/2])
+  vv = np.hstack([v-vd/2.,v[-1]+vd[-1]/2.])
+  if oversample>1:
+    vvo = []
+    for i in range(len(vv)-1):
+      vvo.append(np.linspace(vv[i],vv[i+1],oversample+1)[:-1])
+    vvo.append(vv[-1])
+    vv = np.array(np.hstack(vvo))
+  return vv
+
+def subtractBkg(imgs,nPix=100,dKtype='corners'):
+  """ Opals tend to have different backgroud for every quadrant """
+  if dKtype is 'corners':
+    if imgs.ndim == 2: imgs = imgs[np.newaxis,:]
+    imgs = imgs.astype(np.float)
+    q1 = imgs[:,:nPix,:nPix].mean(-1).mean(-1)
+    imgs[:,:512,:512]-=q1[:,np.newaxis,np.newaxis]
+    q2 = imgs[:,:nPix,-nPix:].mean(-1).mean(-1)
+    imgs[:,:512,-512:]-=q2[:,np.newaxis,np.newaxis]
+    q3 = imgs[:,-nPix:,-nPix:].mean(-1).mean(-1)
+    imgs[:,-512:,-512:]-=q3[:,np.newaxis,np.newaxis]
+    q4 = imgs[:,-nPix:,:nPix].mean(-1).mean(-1)
+    imgs[:,-512:,:512]-=q4[:,np.newaxis,np.newaxis]
+  elif dKtype is 'stripes':
+    if imgs.ndim == 2: imgs = imgs[np.newaxis,:]
+    imgs = imgs.astype(np.float)
+    s1 =  imgs[:,:nPix,:].mean(-2)
+
+
+
+
+  return np.squeeze(imgs)
+
+def getData():
+  t0 = time.time()
+  context = zmq.Context()
+  socket = context.socket(zmq.REQ)
+  socket.connect('tcp://daq-xpp-mon06:12322')
+  #socket.setsockopt(zmq.SUBSCRIBE, b'')
+  while True:
+    socket.send(b"Request")
+    ret    = socket.recv()
+    ret = pickle.loads(ret, encoding='latin1')
+    print('received',ret.keys(),time.time()-t0)
+    t0 = time.time()
+
+if __name__ == "__main__":
+  getData()
+
+class SpecHandler(object):
+  def __init__(self,connectString='tcp://daq-xpp-mon06:12322',spec1name='opal0',spec2name='opal1',roi=[0,1024,0,1024]):
+    self.connectString = connectString
+    self.spec1name = spec1name
+    self.spec2name = spec2name
+    self.dataCollector = []
+    #self.surveyplot = Surveyplot(spec1name=spec1name,spec2name=spec2name)
+    self.roi = roi
+    self.projsimple = ProjSimple(spec1name=spec1name,spec2name=spec2name,dataCollector=self.dataCollector)
+    self._rawContinuousTime = None
+    self.lastDat = None
+    self.openSocket()
+    self.runningPlot = RunningPlot(self.dataCollector)
+
+
+  def openSocket(self):
+    context = zmq.Context()
+    self.context = context
+    self.socket = context.socket(zmq.REQ)
+    self.socket.connect(self.connectString)
+    #self.socket.setsockopt(zmq.SUBSCRIBE, b'')
+
+  def closeSocket(self):
+    del self.context
+    del self.socket
+
+  def getData(self):
+    self.socket.send(b"Request")
+    ret    = self.socket.recv()
+    ret = pickle.loads(ret, encoding='latin1')
+    for sn in [self.spec1name,self.spec2name]:
+      ret[sn] = np.squeeze(alignment.subtractBkg(ret[sn], nPix=100, bkg_type='line'))
+    self.lastDat = ret
+    return ret
+
+  def getRaw(self,repoenConnection=False,doAlign=False,show=False,doFit=False,updateImg=True,updateProj=True,updateFom=True, ispaceCorr=False,flipit=False,xshift=None):
+    if doFit is True: doFit='iminuit'
+
+
+    if repoenConnection:
+      self.closeSocket()
+      self.openSocket()
+    dat = self.getData()
+    im1 = dat[self.spec1name]; im2 = dat[self.spec2name]
+    if doAlign:
+      #t = np.load("gui_align_transform_xppl3716.npy").item()
+      if hasattr(self,'transformer'):
+        algn = self.transformer
+      else:
+        algn = alignment.loadAlignment('last_trafo.npy')
+      t = algn['transform']
+      roi1 = algn['roi1']
+      roi2 = algn['roi2']
+      im1 = im1[roi1]; 
+      if flipit:
+        im1 = im1[::-1]
+      im2 = im2[roi2]
+      alignment.clearCache() 
+      r = alignment.doShot( im1,im2, t, show=show, doFit=doFit, xshift=xshift)
+      self.transformer = dict(transform=r.final_transform,roi1=roi1,roi2=roi2)
+      alignment.saveAlignment('last_trafo.npy',r.final_transform,roi1,roi2)
+
+      im1 = r.im1; im2 = r.im2
+      showDiff = True
+      showRatio = True
+    else:
+      showDiff = False
+      showRatio = False
+    if doAlign:
+      thres = 0.05
+      #im1[im1<thres*np.max(im1.ravel())] = np.nan
+      #im2[im2<thres*np.max(im2.ravel())] = np.nan
+      self.dataCollector.append(\
+	  dict( time=datetime.datetime.now(),
+	    fom=r.fom, 
+	    ratProj = np.nansum(im2/im1,axis=0),
+	    im1Proj = np.nansum(im1,axis=0),
+	    im2Proj = np.nansum(im2,axis=0)))
+      if updateFom:
+        self.runningPlot.updatePlot()
+    #if updateImg:
+      #self.surveyplot.plotImages(im1,im2,showDiff=showDiff,showRatio=showRatio)
+    #if updateProj:
+    self.projsimple.plotProfiles(im1,im2)
+      
+
+  def alignFeatures(self,flipit=False):
+    im1 = copy.copy(self.lastDat[self.spec1name])
+    im2 = copy.copy(self.lastDat[self.spec2name])
+    if flipit:
+      im1 = im1[::-1]
+    roi1 = alignment.findRoi(im1)
+    roi2 = alignment.findRoi(im2)
+    tra = alignment.GuiAlignment(im1[roi1,:],im2[roi2,:],autostart=False)
+    self.transformer = dict(transform=tra.start(),roi1=roi1,roi2=roi2)
+
+  
+  def getRawContinuuous(self,sleepTime,**kwargs):
+    self._rawContinuousTime = sleepTime
+    if not hasattr(self,'_rawContinuousThread'):
+      def upd():
+        while not self._rawContinuousTime is None:
+          self.getRaw(**kwargs)
+          plt.draw()
+          time.sleep(self._rawContinuousTime)
+
+      self._rawContinuousThread = threading.Thread(target=upd)
+      self._rawContinuousThread.start()
+
+
+
+class Surveyplot(object):
+  def __init__(self,spec1name='spec1',spec2name='spec2'):
+    self.fig,self.axs = plt.subplots(4,1,sharex=True,sharey=True)
+    self.axs[0].set_title(spec1name)
+    self.axs[1].set_title(spec2name)
+    self.axs[2].set_title("Difference")
+    self.axs[3].set_title("Ratio")
+
+  def plotImages(self,img1,img2,showDiff=False,showRatio=False):
+    if hasattr(self,'i1'):
+      self.i1.set_data(img1)
+    else:
+      self.i1 = self.axs[0].imshow(img1,origin='lower',interpolation='none')
+    if hasattr(self,'i2'):
+      self.i2.set_data(img2)
+    else:
+      self.i2 = self.axs[1].imshow(img2,origin='lower',interpolation='none')
+    if showDiff:
+      tdiff = img2-img1
+      if hasattr(self,'idiff'):
+        self.idiff.set_data(tdiff)
+      else:
+        self.idiff = self.axs[2].imshow(tdiff,interpolation='none',origin='lower')
+        lms = np.percentile(tdiff,[30,70])
+        self.idiff.set_clim(lms)
+    if showRatio:
+      tratio = img2/img1
+      if hasattr(self,'iratio'):
+        self.iratio.set_data(tratio)
+      else:
+        self.iratio = self.axs[3].imshow(tratio,interpolation='none',origin='lower')
+        lms = np.percentile(tratio,[30,70])
+        self.iratio.set_clim(lms)
+
+class ProjSimple(object):
+  def __init__(self,spec1name='spec1',spec2name='spec2',roi1=[0,1024,0,1024],roi2=[0,1024,0,1024],dataCollector=[]):
+    self.fig,self.axs = plt.subplots(2,1,sharex=True)
+    self.roi1 = roi1
+    self.roi2 = roi2
+    self.spec1name=spec1name
+    self.spec2name=spec2name
+    self.dataCollector = dataCollector
+  
+  #def getROI(self,specNo):
+  def _roiit(self,img,roi):
+    return img[roi[0]:roi[1],roi[2]:roi[3]]
+  
+  def plotProfiles(self,img1,img2):
+    prof1 = np.nansum(self._roiit(img1,self.roi1),0)
+    prof2 = np.nansum(self._roiit(img2,self.roi2),0)
+    if hasattr(self,'l1'):
+      self.l1.set_ydata(prof1)
+    else:
+      self.l1 = self.axs[0].plot(prof1,label=self.spec1name)[0]
+    if hasattr(self,'l2'):
+      self.l2.set_ydata(prof2)
+    else:
+      self.l2 = self.axs[0].plot(prof2,label=self.spec2name)[0]
+      plt.legend()
+
+    if hasattr(self,'lrat'):
+      self.lrat.set_ydata(prof2/prof1)
+    else:
+      self.lrat = self.axs[1].plot(prof2/prof1,'k',label='ratio')[0]
+      self.axs[1].set_ylim(0,2)
+    if len(self.dataCollector) > 0 :
+      im1Proj = np.asarray([i['im1Proj'] for i in self.dataCollector])
+      im2Proj = np.asarray([i['im2Proj'] for i in self.dataCollector])
+      #print(im1Proj.shape,im2Proj.shape)
+      ratAv = np.median(im2Proj[-10:],0)/np.median(im1Proj[-10:],0)
+      if hasattr(self,'lratAv'):
+        self.lratAv.set_ydata(ratAv)
+      else:
+        self.lratAv = self.axs[1].plot(ratAv,'r',label='ratio Avg')[0]
+        self.axs[1].set_ylim(0,2)
+
+
+class RunningPlot(object):
+  def __init__(self,dataCollector):
+    self.dataCollector = dataCollector
+    self.fig,self.axs = plt.subplots(1,1)
+    
+  def updatePlot(self):
+    if len(self.dataCollector)>0:
+      times = np.asarray(([i['time'] for i in self.dataCollector]))
+      foms = np.asarray(([i['fom'] for i in self.dataCollector]))
+      im1Proj = np.asarray(([i['im1Proj'] for i in self.dataCollector]))
+      im2Proj = np.asarray(([i['im2Proj'] for i in self.dataCollector]))
+
+      if hasattr(self,'fomline'):
+        self.fomline.set_ydata(foms)
+        self.fomline.set_xdata(times)
+        self.axs.set_xlim(np.min(times)-datetime.timedelta(0,10),np.max(times)+datetime.timedelta(0,10))
+
+	#self.axs[0].autoscale(enable=True,axis='x')
+
+      else:
+        self.fomline = self.axs.plot(times,foms,'o-')[0]
+        self.axs.autoscale(enable=True,axis='x')
+      #if hasattr(self,'ratioimg'):
+	#self.ratioimg.set_data()
+	#self.ratioimg.set_ydata(foms)
+      #else:
+	#self.axs[0].plot(times,foms,'o-')
+
+
+   
+
+  
+  #def plotOrUpdate(img1,img2):
+    #if hasattr(self,i1):
+      #self.i1.set_data(img1)
+    #else:
+      #self.i1 = self.axs.imshow(img1,interpolate='none',origin='bottom')
+    #if hasattr(self,i1):
+      #self.i1.set_data(img1)
+    #else:
+      #self.i1 = self.axs.imshow(img1,interpolate='none',origin='bottom')
+

online/sMhelper_dev_pg.py

@@ -0,0 +1,345 @@
+import sys
+import zmq
+import numpy as np
+import time
+import pickle
+import alignment
+import matplotlib.pyplot as plt
+import threading
+import datetime
+import copy
+import pyqtgraph as pg
+
+plt.rcParams['image.cmap'] = 'viridis'
+
+def histVec(v,oversample=1):
+  v = np.atleast_1d(v)
+  v = np.unique(v)
+  vd = np.diff(v)
+  vd = np.hstack([vd[0],vd])
+  #vv = np.hstack([v-vd/2,v[-1]+vd[-1]/2])
+  vv = np.hstack([v-vd/2.,v[-1]+vd[-1]/2.])
+  if oversample>1:
+    vvo = []
+    for i in range(len(vv)-1):
+      vvo.append(np.linspace(vv[i],vv[i+1],oversample+1)[:-1])
+    vvo.append(vv[-1])
+    vv = np.array(np.hstack(vvo))
+  return vv
+
+def subtractBkg(imgs,nPix=100,dKtype='corners'):
+  """ Opals tend to have different backgroud for every quadrant """
+  if dKtype is 'corners':
+    if imgs.ndim == 2: imgs = imgs[np.newaxis,:]
+    imgs = imgs.astype(np.float)
+    q1 = imgs[:,:nPix,:nPix].mean(-1).mean(-1)
+    imgs[:,:512,:512]-=q1[:,np.newaxis,np.newaxis]
+    q2 = imgs[:,:nPix,-nPix:].mean(-1).mean(-1)
+    imgs[:,:512,-512:]-=q2[:,np.newaxis,np.newaxis]
+    q3 = imgs[:,-nPix:,-nPix:].mean(-1).mean(-1)
+    imgs[:,-512:,-512:]-=q3[:,np.newaxis,np.newaxis]
+    q4 = imgs[:,-nPix:,:nPix].mean(-1).mean(-1)
+    imgs[:,-512:,:512]-=q4[:,np.newaxis,np.newaxis]
+  elif dKtype is 'stripes':
+    if imgs.ndim == 2: imgs = imgs[np.newaxis,:]
+    imgs = imgs.astype(np.float)
+    s1 =  imgs[:,:nPix,:].mean(-2)
+
+
+
+
+  return np.squeeze(imgs)
+
+def getData():
+  t0 = time.time()
+  context = zmq.Context()
+  socket = context.socket(zmq.REQ)
+  socket.connect('tcp://daq-xpp-mon06:12322')
+  #socket.setsockopt(zmq.SUBSCRIBE, b'')
+  while True:
+    socket.send(b"Request")
+    ret    = socket.recv()
+    ret = pickle.loads(ret, encoding='latin1')
+    print('received',ret.keys(),time.time()-t0)
+    t0 = time.time()
+
+if __name__ == "__main__":
+  getData()
+
+class SpecHandler(object):
+  def __init__(self,connectString='tcp://daq-xpp-mon06:12322',spec1name='opal0',spec2name='opal1',roi=[0,1024,0,1024]):
+    self.connectString = connectString
+    self.spec1name = spec1name
+    self.spec2name = spec2name
+    self.dataCollector = []
+    #self.surveyplot = Surveyplot(spec1name=spec1name,spec2name=spec2name)
+    self.roi = roi
+    self.projsimple = ProjSimple(spec1name=spec1name,spec2name=spec2name,dataCollector=self.dataCollector)
+    self._rawContinuousTime = None
+    self.lastDat = None
+    self.openSocket()
+    #self.runningPlot = RunningPlot(self.dataCollector)
+
+
+  def openSocket(self):
+    context = zmq.Context()
+    self.context = context
+    self.socket = context.socket(zmq.REQ)
+    self.socket.connect(self.connectString)
+    #self.socket.setsockopt(zmq.SUBSCRIBE, b'')
+
+  def closeSocket(self):
+    del self.context
+    del self.socket
+
+  def getData(self):
+    self.socket.send(b"Request")
+    ret    = self.socket.recv()
+    ret = pickle.loads(ret, encoding='latin1')
+    for sn in [self.spec1name,self.spec2name]:
+      ret[sn] = np.squeeze(alignment.subtractBkg(ret[sn], nPix=100, bkg_type='line'))
+    self.lastDat = ret
+    return ret
+
+  def getRaw(self,repoenConnection=False,doAlign=False,show=False,doFit=False,updateImg=True,updateProj=True,updateFom=True):
+    if doFit is True: doFit='iminuit'
+
+
+    if repoenConnection:
+      self.closeSocket()
+      self.openSocket()
+    dat = self.getData()
+    im1 = dat[self.spec1name]; im2 = dat[self.spec2name]
+    if doAlign:
+      #t = np.load("gui_align_transform_xppl3716.npy").item()
+      if hasattr(self,'transformer'):
+        algn = self.transformer
+      else:
+        algn = alignment.loadAlignment('last_trafo.npy')
+      t = algn['transform']
+      roi1 = algn['roi1']
+      roi2 = algn['roi2']
+      
+      r = alignment.doShot( im1[roi1,:],im2[roi2,:],t, show=show, doFit=doFit)
+      self.transformer = dict(transform=r.final_transform,roi1=roi1,roi2=roi2)
+      alignment.saveAlignment('last_trafo.npy',r.final_transform,roi1,roi2)
+
+      im1 = r.im1; im2 = r.im2
+      showDiff = True
+      showRatio = True
+    else:
+      showDiff = False
+      showRatio = False
+    if doAlign:
+      thres = 0.05
+      #im1[im1<thres*np.max(im1.ravel())] = np.nan
+      #im2[im2<thres*np.max(im2.ravel())] = np.nan
+      self.dataCollector.append(\
+	  dict( time=datetime.datetime.now(),
+	    fom=r.fom, 
+	    ratProj = np.nansum(im2/im1,axis=0),
+	    im1Proj = np.nansum(im1,axis=0),
+	    im2Proj = np.nansum(im2,axis=0)))
+	  #if updateFom:
+	#self.runningPlot.updatePlot()
+    #if updateImg:
+      #self.surveyplot.plotImages(im1,im2,showDiff=showDiff,showRatio=showRatio)
+    #if updateProj:
+    self.projsimple.plotProfiles(im1,im2)
+      
+
+  def alignFeatures(self):
+    im1 = copy.copy(self.lastDat[self.spec1name])
+    im2 = copy.copy(self.lastDat[self.spec2name])
+    roi1 = alignment.findRoi(im1)
+    roi2 = alignment.findRoi(im2)
+
+    tra = alignment.GuiAlignment(im1[roi1,:],im2[roi2,:],autostart=False)
+    self.transformer = dict(transform=tra.start(),roi1=roi1,roi2=roi2)
+
+  
+  def getRawContinuuous(self,sleepTime,**kwargs):
+    self._rawContinuousTime = sleepTime
+    if not hasattr(self,'_rawContinuousThread'):
+      def upd():
+        while not self._rawContinuousTime is None:
+          self.getRaw(**kwargs)
+          plt.draw()
+          time.sleep(self._rawContinuousTime)
+
+      self._rawContinuousThread = threading.Thread(target=upd)
+      self._rawContinuousThread.start()
+
+
+
+class Surveyplot(object):
+  def __init __(self,spec1name='spec1',spec2name='spec2'):
+    self.fig,self.axs = plt.subplots(4,1,sharex=True,sharey=True)
+    self.axs[0].set_title(spec1name)
+    self.axs[1].set_title(spec2name)
+    self.axs[2].set_title("Difference")
+    self.axs[3].set_title("Ratio")
+
+  def plotImages(self,img1,img2,showDiff=False,showRatio=False):
+    if hasattr(self,'i1'):
+      self.i1.set_data(img1)
+    else:
+      self.i1 = self.axs[0].imshow(img1,origin='lower',interpolation='none')
+    if hasattr(self,'i2'):
+      self.i2.set_data(img2)
+    else:
+      self.i2 = self.axs[1].imshow(img2,origin='lower',interpolation='none')
+    if showDiff:
+      tdiff = img2-img1
+      if hasattr(self,'idiff'):
+        self.idiff.set_data(tdiff)
+      else:
+        self.idiff = self.axs[2].imshow(tdiff,interpolation='none',origin='lower')
+        lms = np.percentile(tdiff,[30,70])
+        self.idiff.set_clim(lms)
+    if showRatio:
+      tratio = img2/img1
+      if hasattr(self,'iratio'):
+        self.iratio.set_data(tratio)
+      else:
+        self.iratio = self.axs[3].imshow(tratio,interpolation='none',origin='lower')
+        lms = np.percentile(tratio,[30,70])
+        self.iratio.set_clim(lms)
+
+class ProjSimple(object):
+  def __init__(self,spec1name='sp ec1',spec2name='spec2',roi1=[0,1024,0,1024],roi2=[0,1024,0,1024],dataCollector=[]):
+    self.fig,self.axs = plt.subplots(2,1,sharex=True)
+    self.roi1 = roi1
+    self.roi2 = roi2
+    self.spec1name=spec1name
+    self.spec2name=spec2name
+    self.dataCollector = dataCollector
+  
+  #def getROI(self,specNo):
+  def _roiit(self,img,roi):
+    return img[roi[0]:roi[1],roi[2]:roi[3]]
+  
+  def plotProfiles(self,img1,img2):
+    prof1 = np.nansum(self._roiit(img1,self.roi1),0)
+    prof2 = np.nansum(self._roiit(img2,self.roi2),0)
+    if hasattr(self,'l1'):
+      self.l1.set_ydata(prof1)
+    else:
+      self.l1 = self.axs[0].plot(prof1,label=self.spec1name)[0]
+    if hasattr(self,'l2'):
+      self.l2.set_ydata(prof2)
+    else:
+      self.l2 = self.axs[0].plot(prof2,label=self.spec2name)[0]
+      plt.legend()
+
+    if hasattr(self,'lrat'):
+      self.lrat.set_ydata(prof2/prof1)
+    else:
+      self.lrat = self.axs[1].plot(prof2/prof1,'k',label='ratio')[0]
+      self.axs[1].set_ylim(0,2)
+    if len(self.dataCollector) > 0 :
+      im1Proj = np.asarray([i['im1Proj'] for i in self.dataCollector])
+      im2Proj = np.asarray([i['im2Proj'] for i in self.dataCollector])
+      #print(im1Proj.shape,im2Proj.shape)
+      ratAv = np.median(im2Proj,0)/np.median(im1Proj,0)
+      if hasattr(self,'lratAv'):
+        self.lratAv.set_ydata(ratAv)
+      else:
+        self.lratAv = self.axs[1].plot(ratAv,'r',label='ratio Avg')[0]
+        self.axs[1].set_ylim(0,2)
+
+class ProjSimplePG(object):
+  def __init__(self,spec1name='sp ec1',spec2name='spec2',roi1=[0,1024,0,1024],roi2=[0,1024,0,1024],dataCollector=[]):
+    from pyqtgraph.Qt import QtGui, QtCore
+    import numpy as np
+    import pyqtgraph as pg
+    self.app = QtGui.Qapplication([])
+    self.win = pg.GraphicsWindow(title='Single shot XANES')
+    pg.setConfigOptions(antialias=True)
+    self.sp1 = win.addPlot(title='spectrometer spectra')
+    self.sp2 = win.addPlot(title='Ratio')
+
+
+    self.fig,self.axs = plt.subplots(2,1,sharex=True)
+
+    self.roi1 = roi1
+    self.roi2 = roi2
+    self.spec1name=spec1name
+    self.spec2name=spec2name
+    self.dataCollector = dataCollector
+  
+  #def getROI(self,specNo):
+  def _roiit(self,img,roi):
+    return img[roi[0]:roi[1],roi[2]:roi[3]]
+  
+  def plotProfiles(self,img1,img2):
+    prof1 = np.nansum(self._roiit(img1,self.roi1),0)
+    prof2 = np.nansum(self._roiit(img2,self.roi2),0)
+    if hasattr(self,'l1'):
+      self.l1.set_ydata(prof1)
+    else:
+      self.l1 = self.axs[0].plot(prof1,label=self.spec1name)[0]
+    if hasattr(self,'l2'):
+      self.l2.set_ydata(prof2)
+    else:
+      self.l2 = self.axs[0].plot(prof2,label=self.spec2name)[0]
+      plt.legend()
+
+    if hasattr(self,'lrat'):
+      self.lrat.set_ydata(prof2/prof1)
+    else:
+      self.lrat = self.axs[1].plot(prof2/prof1,'k',label='ratio')[0]
+      self.axs[1].set_ylim(0,2)
+    if len(self.dataCollector) > 0 :
+      im1Proj = np.asarray([i['im1Proj'] for i in self.dataCollector])
+      im2Proj = np.asarray([i['im2Proj'] for i in self.dataCollector])
+      #print(im1Proj.shape,im2Proj.shape)
+      ratAv = np.median(im2Proj,0)/np.median(im1Proj,0)
+      if hasattr(self,'lratAv'):
+        self.lratAv.set_ydata(ratAv)
+      else:
+        self.lratAv = self.axs[1].plot(ratAv,'r',label='ratio Avg')[0]
+        self.axs[1].set_ylim(0,2)
+
+class RunningPlot(object):
+  def __init__(self,dataCollector):
+    self.dataCollector = dataCollector
+    self.fig,self.axs = plt.subplots(1,1)
+    
+  def updatePlot(self):
+    if len(self.dataCollector)>0:
+      times = np.asarray(([i['time'] for i in self.dataCollector]))
+      foms = np.asarray(([i['fom'] for i in self.dataCollector]))
+      im1Proj = np.asarray(([i['im1Proj'] for i in self.dataCollector]))
+      im2Proj = np.asarray(([i['im2Proj'] for i in self.dataCollector]))
+
+      if hasattr(self,'fomline'):
+        self.fomline.set_ydata(foms)
+        self.fomline.set_xdata(times)
+        self.axs.set_xlim(np.min(times)-datetime.timedelta(0,10),np.max(times)+datetime.timedelta(0,10))
+
+	#self.axs[0].autoscale(enable=True,axis='x')
+
+      else:
+        self.fomline = self.axs.plot(times,foms,'o-')[0]
+        self.axs.autoscale(enable=True,axis='x')
+      #if hasattr(self,'ratioimg'):
+	#self.ratioimg.set_data()
+	#self.ratioimg.set_ydata(foms)
+      #else:
+	#self.axs[0].plot(times,foms,'o-')
+
+
+   
+
+  
+  #def plotOrUpdate(img1,img2):
+    #if hasattr(self,i1):
+      #self.i1.set_data(img1)
+    #else:
+      #self.i1 = self.axs.imshow(img1,interpolate='none',origin='bottom')
+    #if hasattr(self,i1):
+      #self.i1.set_data(img1)
+    #else:
+      #self.i1 = self.axs.imshow(img1,interpolate='none',origin='bottom')
+

online/sMhelper_upd.py

@@ -0,0 +1,296 @@
+import sys
+import zmq
+import numpy as np
+import time
+import pickle
+import alignment
+import matplotlib.pyplot as plt
+import threading
+import datetime
+import copy
+
+plt.rcParams['image.cmap'] = 'viridis'
+
+def histVec(v,oversample=1):
+  v = np.atleast_1d(v)
+  v = np.unique(v)
+  vd = np.diff(v)
+  vd = np.hstack([vd[0],vd])
+  #vv = np.hstack([v-vd/2,v[-1]+vd[-1]/2])
+  vv = np.hstack([v-vd/2.,v[-1]+vd[-1]/2.])
+  if oversample>1:
+    vvo = []
+    for i in range(len(vv)-1):
+      vvo.append(np.linspace(vv[i],vv[i+1],oversample+1)[:-1])
+    vvo.append(vv[-1])
+    vv = np.array(np.hstack(vvo))
+  return vv
+
+def subtractBkg(imgs,nPix=100,dKtype='corners'):
+  """ Opals tend to have different backgroud for every quadrant """
+  if dKtype is 'corners':
+    if imgs.ndim == 2: imgs = imgs[np.newaxis,:]
+    imgs = imgs.astype(np.float)
+    q1 = imgs[:,:nPix,:nPix].mean(-1).mean(-1)
+    imgs[:,:512,:512]-=q1[:,np.newaxis,np.newaxis]
+    q2 = imgs[:,:nPix,-nPix:].mean(-1).mean(-1)
+    imgs[:,:512,-512:]-=q2[:,np.newaxis,np.newaxis]
+    q3 = imgs[:,-nPix:,-nPix:].mean(-1).mean(-1)
+    imgs[:,-512:,-512:]-=q3[:,np.newaxis,np.newaxis]
+    q4 = imgs[:,-nPix:,:nPix].mean(-1).mean(-1)
+    imgs[:,-512:,:512]-=q4[:,np.newaxis,np.newaxis]
+  elif dKtype is 'stripes':
+    if imgs.ndim == 2: imgs = imgs[np.newaxis,:]
+    imgs = imgs.astype(np.float)
+    s1 =  imgs[:,:nPix,:].mean(-2)
+
+
+
+
+  return np.squeeze(imgs)
+
+def getData():
+  t0 = time.time()
+  context = zmq.Context()
+  socket = context.socket(zmq.REQ)
+  socket.connect('tcp://daq-xpp-mon06:12322')
+  #socket.setsockopt(zmq.SUBSCRIBE, b'')
+  while True:
+    socket.send(b"Request")
+    ret    = socket.recv()
+    ret = pickle.loads(ret, encoding='latin1')
+    print('received',ret.keys(),time.time()-t0)
+    t0 = time.time()
+
+if __name__ == "__main__":
+  getData()
+
+class SpecHandler(object):
+  def __init__(self,connectString='tcp://daq-xpp-mon06:12322',spec1name='opal0',spec2name='opal1',roi=[0,1024,0,1024]):
+    self.connectString = connectString
+    self.spec1name = spec1name
+    self.spec2name = spec2name
+    self.dataCollector = []
+    #self.surveyplot = Surveyplot(spec1name=spec1name,spec2name=spec2name)
+    self.roi = roi
+    self.projsimple = ProjSimple(spec1name=spec1name,spec2name=spec2name,dataCollector=self.dataCollector)
+    self._rawContinuousTime = None
+    self.lastDat = None
+    self.openSocket()
+    #self.runningPlot = RunningPlot(self.dataCollector)
+
+
+  def openSocket(self):
+    context = zmq.Context()
+    self.context = context
+    self.socket = context.socket(zmq.REQ)
+    self.socket.connect(self.connectString)
+    #self.socket.setsockopt(zmq.SUBSCRIBE, b'')
+
+  def closeSocket(self):
+    del self.context
+    del self.socket
+
+  def getData(self):
+    self.socket.send(b"Request")
+    ret    = self.socket.recv()
+    ret = pickle.loads(ret, encoding='latin1')
+    for sn in [self.spec1name,self.spec2name]:
+      ret[sn] = np.squeeze(alignment.subtractBkg(ret[sn], nPix=100, bkg_type='line'))
+    self.lastDat = ret
+    return ret
+
+  def getRaw(self,repoenConnection=False,doAlign=False,show=False,doFit=False,updateImg=True,updateProj=True,updateFom=True,flipit=False):
+    if doFit is True: doFit='iminuit'
+
+
+    if repoenConnection:
+      self.closeSocket()
+      self.openSocket()
+    dat = self.getData()
+    im1 = dat[self.spec1name]; im2 = dat[self.spec2name]
+    if doAlign:
+      #t = np.load("gui_align_transform_xppl3716.npy").item()
+      #if hasattr(self,'transformer'):
+	#algn = self.transformer
+      #else:
+      algn = alignment.loadAlignment('last_trafo.npy')
+      t = algn['transform']
+      roi1 = algn['roi1']
+      roi2 = algn['roi2']
+      im1 = im1[roi1]; 
+      if flipit:
+	im1 = im1[::-1]
+      im2 = im2[roi2]
+      
+      r = alignment.doShot( im1,im2, t, show=show, doFit=doFit)
+      self.transformer = dict(transform=r.final_transform,roi1=roi1,roi2=roi2)
+      alignment.saveAlignment('last_trafo.npy',r.final_transform,roi1,roi2)
+
+      im1 = r.im1; im2 = r.im2
+      showDiff = True
+      showRatio = True
+    else:
+      showDiff = False
+      showRatio = False
+    if doAlign:
+      thres = 0.05
+      #im1[im1<thres*np.max(im1.ravel())] = np.nan
+      #im2[im2<thres*np.max(im2.ravel())] = np.nan
+      self.dataCollector.append(\
+	  dict( time=datetime.datetime.now(),
+	    fom=r.fom, 
+	    ratProj = np.nansum(im2/im1,axis=0),
+	    im1Proj = np.nansum(im1,axis=0),
+	    im2Proj = np.nansum(im2,axis=0)))
+	  #if updateFom:
+	#self.runningPlot.updatePlot()
+    #if updateImg:
+      #self.surveyplot.plotImages(im1,im2,showDiff=showDiff,showRatio=showRatio)
+    #if updateProj:
+    self.projsimple.plotProfiles(im1,im2)
+      
+
+  def alignFeatures(self):
+    im1 = copy.copy(self.lastDat[self.spec1name])
+    im2 = copy.copy(self.lastDat[self.spec2name])
+    roi1 = alignment.findRoi(im1)
+    roi2 = alignment.findRoi(im2)
+
+    tra = alignment.GuiAlignment(im1[roi1,:],im2[roi2,:],autostart=False)
+    self.transformer = dict(transform=tra.start(),roi1=roi1,roi2=roi2)
+
+  
+  def getRawContinuuous(self,sleepTime,**kwargs):
+    self._rawContinuousTime = sleepTime
+    if not hasattr(self,'_rawContinuousThread'):
+      def upd():
+        while not self._rawContinuousTime is None:
+          self.getRaw(**kwargs)
+          plt.draw()
+          time.sleep(self._rawContinuousTime)
+
+      self._rawContinuousThread = threading.Thread(target=upd)
+      self._rawContinuousThread.start()
+
+
+
+class Surveyplot(object):
+  def __init__(self,spec1name='spec1',spec2name='spec2'):
+    self.fig,self.axs = plt.subplots(4,1,sharex=True,sharey=True)
+    self.axs[0].set_title(spec1name)
+    self.axs[1].set_title(spec2name)
+    self.axs[2].set_title("Difference")
+    self.axs[3].set_title("Ratio")
+
+  def plotImages(self,img1,img2,showDiff=False,showRatio=False):
+    if hasattr(self,'i1'):
+      self.i1.set_data(img1)
+    else:
+      self.i1 = self.axs[0].imshow(img1,origin='lower',interpolation='none')
+    if hasattr(self,'i2'):
+      self.i2.set_data(img2)
+    else:
+      self.i2 = self.axs[1].imshow(img2,origin='lower',interpolation='none')
+    if showDiff:
+      tdiff = img2-img1
+      if hasattr(self,'idiff'):
+        self.idiff.set_data(tdiff)
+      else:
+        self.idiff = self.axs[2].imshow(tdiff,interpolation='none',origin='lower')
+        lms = np.percentile(tdiff,[30,70])
+        self.idiff.set_clim(lms)
+    if showRatio:
+      tratio = img2/img1
+      if hasattr(self,'iratio'):
+        self.iratio.set_data(tratio)
+      else:
+        self.iratio = self.axs[3].imshow(tratio,interpolation='none',origin='lower')
+        lms = np.percentile(tratio,[30,70])
+        self.iratio.set_clim(lms)
+
+class ProjSimple(object):
+  def __init__(self,spec1name='spec1',spec2name='spec2',roi1=[0,1024,0,1024],roi2=[0,1024,0,1024],dataCollector=[]):
+    self.fig,self.axs = plt.subplots(2,1,sharex=True)
+    self.roi1 = roi1
+    self.roi2 = roi2
+    self.spec1name=spec1name
+    self.spec2name=spec2name
+    self.dataCollector = dataCollector
+  
+  #def getROI(self,specNo):
+  def _roiit(self,img,roi):
+    return img[roi[0]:roi[1],roi[2]:roi[3]]
+  
+  def plotProfiles(self,img1,img2):
+    prof1 = np.nansum(self._roiit(img1,self.roi1),0)
+    prof2 = np.nansum(self._roiit(img2,self.roi2),0)
+    if hasattr(self,'l1'):
+      self.l1.set_ydata(prof1)
+    else:
+      self.l1 = self.axs[0].plot(prof1,label=self.spec1name)[0]
+    if hasattr(self,'l2'):
+      self.l2.set_ydata(prof2)
+    else:
+      self.l2 = self.axs[0].plot(prof2,label=self.spec2name)[0]
+      plt.legend()
+
+    if hasattr(self,'lrat'):
+      self.lrat.set_ydata(prof2/prof1)
+    else:
+      self.lrat = self.axs[1].plot(prof2/prof1,'k',label='ratio')[0]
+      self.axs[1].set_ylim(0,2)
+    if len(self.dataCollector) > 0 :
+      im1Proj = np.asarray([i['im1Proj'] for i in self.dataCollector])
+      im2Proj = np.asarray([i['im2Proj'] for i in self.dataCollector])
+      #print(im1Proj.shape,im2Proj.shape)
+      ratAv = np.median(im2Proj[-10:,:],0)/np.median(im1Proj[-10:,:],0)
+      if hasattr(self,'lratAv'):
+        self.lratAv.set_ydata(ratAv)
+      else:
+        self.lratAv = self.axs[1].plot(ratAv,'r',label='ratio Avg')[0]
+        self.axs[1].set_ylim(0,2)
+
+
+class RunningPlot(object):
+  def __init__(self,dataCollector):
+    self.dataCollector = dataCollector
+    self.fig,self.axs = plt.subplots(1,1)
+    
+  def updatePlot(self):
+    if len(self.dataCollector)>0:
+      times = np.asarray(([i['time'] for i in self.dataCollector]))
+      foms = np.asarray(([i['fom'] for i in self.dataCollector]))
+      im1Proj = np.asarray(([i['im1Proj'] for i in self.dataCollector]))
+      im2Proj = np.asarray(([i['im2Proj'] for i in self.dataCollector]))
+
+      if hasattr(self,'fomline'):
+        self.fomline.set_ydata(foms)
+        self.fomline.set_xdata(times)
+        self.axs.set_xlim(np.min(times)-datetime.timedelta(0,10),np.max(times)+datetime.timedelta(0,10))
+
+	#self.axs[0].autoscale(enable=True,axis='x')
+
+      else:
+        self.fomline = self.axs.plot(times,foms,'o-')[0]
+        self.axs.autoscale(enable=True,axis='x')
+      #if hasattr(self,'ratioimg'):
+	#self.ratioimg.set_data()
+	#self.ratioimg.set_ydata(foms)
+      #else:
+	#self.axs[0].plot(times,foms,'o-')
+
+
+   
+
+  
+  #def plotOrUpdate(img1,img2):
+    #if hasattr(self,i1):
+      #self.i1.set_data(img1)
+    #else:
+      #self.i1 = self.axs.imshow(img1,interpolate='none',origin='bottom')
+    #if hasattr(self,i1):
+      #self.i1.set_data(img1)
+    #else:
+      #self.i1 = self.axs.imshow(img1,interpolate='none',origin='bottom')
+

online/server.py

@@ -0,0 +1,41 @@
+import zmq
+import time
+import psana
+import numpy as np
+
+context = zmq.Context()
+socket = context.socket(zmq.REP)
+socket.bind("tcp://*:12322")
+poller = zmq.Poller()
+poller.register(socket, zmq.POLLIN)
+
+#run = 8
+#ds = psana.DataSource('exp=xppl3816:run=%d:smd' %run)
+ds = psana.DataSource('shmem=XPP.0:stop=no')
+epics = ds.env().epicsStore()
+opal_0_detector = psana.Detector('opal_0')
+opal_1_detector = psana.Detector('opal_1')
+#ipm3_src = psana.Source('BldInfo(XppSb3_Ipm)')
+
+t0 = time.time()
+for i, evt in enumerate(ds.events()):
+    events = dict(poller.poll(0))
+    if socket in events and events[socket] == zmq.POLLIN:
+        opal_0 = opal_0_detector.raw(evt)
+        opal_1 = opal_1_detector.raw(evt)
+        cntr = ds.env().configStore().get(psana.ControlData.ConfigV3, psana.Source()).pvControls()
+        if len(cntr) > 0:
+          cpvName = cntr[0].name()
+          cpvValue = cntr[0].value()
+        else:
+          cpvName = None
+          cpvValue = None
+    
+        if opal_0 is None or opal_1 is None:
+            print 'missing opal'
+            continue
+        message = socket.recv()
+        socket.send_pyobj( dict( opal0 = opal_0, opal1 = opal_1, cPv = dict(name=cpvName,value=cpvValue)) )
+        print 'Shot',i, 'sent; time since starting:', time.time()-t0
+
+

online/server_pub.py

@@ -0,0 +1,41 @@
+import zmq
+import time
+import psana
+import numpy as np
+
+context = zmq.Context()
+socket = context.socket(zmq.PUB)
+socket.bind("tcp://*:12322")
+
+#run = 8
+#ds = psana.DataSource('exp=xppl3816:run=%d:smd' %run)
+ds = psana.DataSource('shmem=XPP.0:stop=no')
+epics = ds.env().epicsStore()
+opal_0_detector = psana.Detector('opal_0')
+opal_1_detector = psana.Detector('opal_1')
+opal_1_detector = psana.Detector('opal_1')
+#ipm3_src = psana.Source('BldInfo(XppSb3_Ipm)')
+
+t0 = time.time()
+for i, evt in enumerate(ds.events()):
+    #if i % 20 != 0:
+    #    continue
+    opal_0 = opal_0_detector.raw(evt)
+#    opal_2 = np.random.random((1024, 1024))#opal_2_detector.raw(evt)
+    opal_1 = opal_1_detector.raw(evt)
+    cntr = ds.env().configStore().get(psana.ControlData.ConfigV3, psana.Source()).pvControls()
+    if len(cntr) > 0:
+      cpvName = cntr.pvControlsu()[0].name()
+      cpvValue = cntr.pvControls()[0].value()
+    else:
+      cpvName = None
+      cpvValue = None
+    
+    if opal_0 is None or opal_1 is None:
+        print 'missing opal'
+        continue
+    socket.send_pyobj( dict( opal0 = opal_0, opal1 = opal_1, cPv = dict(name=cpvName,value=cpvValue)) )
+    print 'Shot',i, 'sent; time since starting:', time.time()-t0
+    time.sleep(1.0)
+
+

x3py_config

@@ -0,0 +1 @@
+#cachepath = "/reg/d/psdm/xpp/xppl3716/scratch/mc/cache"

xanes_analyzeRun.py

@@ -0,0 +1,291 @@
+import os
+import sys
+import numpy as np
+np.warnings.simplefilter('ignore')
+import time
+import matplotlib.pyplot as plt
+import h5py
+import re
+
+from x3py import x3py
+import alignment
+import mcutils as mc
+
+kw_2dplot = dict( 
+  interpolation = "none",
+  aspect        = "auto",
+  cmap          = plt.cm.viridis
+  )
+
+
+g_exp   = "mecl3616"
+g_exp   = "xppl3716"
+g_bml   = g_exp[:3]
+
+x3py.config.updateBeamline(g_bml)
+
+g_folder_init = g_exp+"_init_pars/"
+g_folder_out  = g_exp+"_output/"
+g_folder_data = "/reg/d/psdm/"+g_bml+"/"+ g_exp +"/hdf5/"
+import socket
+hostname = socket.gethostname()
+if hostname == "x1":
+  g_folder_data = "/home/marco/temp"
+
+# set defaults based on experiment
+if g_bml == "xpp":
+  g_roi_height = 200
+  g_swapx = False
+  g_swapy = False
+else:
+  g_roi_height = 100
+  g_swapx = True
+  g_swapy = False
+
+print("Working on experiment",g_exp,"(beamline %s)"%g_bml)
+print(" folder data      →",g_folder_data)
+print(" folder init_pars →",g_folder_init)
+print(" folder outout    →",g_folder_out)
+
+#g_folder = "/reg/d/psdm/xpp/xppl3716/ftc/hdf5/"
+
+def readDataset(fnameOrRun=7,
+  force=False,
+  doBkgSub=False):
+  if isinstance(fnameOrRun,str) and (fnameOrRun[-3:]=="npz"):
+    d = x3py.toolsVarious.DropObject()
+    temp = np.load(fnameOrRun)
+    spec1 = temp["spec1"]
+    spec2 = temp["spec2"]
+    nS    = spec1.shape[0]
+    d.spec1 = x3py.toolsDetectors.wrapArray("spec1",spec1,time=np.arange(nS))
+    d.spec2 = x3py.toolsDetectors.wrapArray("spec2",spec2,time=np.arange(nS))
+  else:
+    if isinstance(fnameOrRun,int):
+      fnameOrRun=g_folder_data+"/"+g_exp+"-r%04d.h5" % fnameOrRun
+    d = x3py.Dataset(fnameOrRun,detectors=["opal0","opal1","fee_spec","opal2"])
+    if g_bml == "xpp":
+      d.spec1 = d.opal0
+      d.spec2 = d.opal1
+    else:
+      d.spec1 = d.fee_spec
+      d.spec2 = d.opal2
+  if not hasattr(d,"scan"):
+    d.scan = x3py.toolsVarious.DropObject()
+    d.scan.scanmotor0_values = [0,]
+  return d
+
+
+def getCenter(img,axis=0,threshold=0.05):
+  img = img.copy()
+  img[img<img.max()*threshold] = 0
+  if axis == 1: img=img.T
+  p   = img.mean(1)
+  x   = np.arange(img.shape[0])
+  return int(np.sum(x*p)/np.sum(p))
+
+def showShots(im1,im2):
+  nS = im1.shape[0]
+  fig,ax = plt.subplots(2,nS,sharex=True,sharey=True)
+  if im1.ndim == 3:
+    for a,i1,i2 in zip(ax.T,im1,im2):
+      a[0].imshow(i1.T,**kw_2dplot)
+      a[1].imshow(i2.T,**kw_2dplot)
+  else:
+    for a,p1,p2 in zip(ax.T,im1,im2):
+      a[0].plot(p1)
+      a[1].plot(p2)
+  
+class AnalyzeRun(object):
+  def __init__(self,run,initAlign="auto",swapx=g_swapx,swapy=g_swapy):
+    """ swapx → swap x axis of first spectrometer
+        swapy → swap y axis of first spectrometer
+    """
+    self.d = readDataset(run)
+    if isinstance(run,str):
+      run = int( re.search("\d{3,4}",run).group() )
+    self.run = run
+    self.results = dict()
+    self.swap = (swapx,swapy)
+    #self.clearCache()
+
+    d = self.d
+    self.spec1 = d.spec1 ; # spec1 is the one that is moved
+    self.spec2 = d.spec2 ; 
+
+    try:
+      self.loadTransform(initAlign)
+    except (AttributeError,FileNotFoundError):
+      if initAlign is None:
+        print("Set to default transform")
+        self.initAlign = self.setDefaultTransform()
+    #self.initAlign  = initAlign
+
+  def getShot(self,shot=0,calib=None,bkgSub="line",roi=g_roi_height):
+    # read data
+    im1 = self.spec1.getShots(shot,calib=calib)
+    im2 = self.spec2.getShots(shot,calib=calib)
+    # subtractBkg bkg
+    im1 = alignment.subtractBkg(im1,bkg_type=bkgSub)
+    im2 = alignment.subtractBkg(im2,bkg_type=bkgSub)
+    # rebin and swap im1 if necessary
+    if im1.shape[-1] != 1024:
+      im1 = mc.rebin(im1, (im1.shape[0],im1.shape[1],1024) )
+    if self.swap[0]:
+      im1 = im1[:,:,::-1]
+    if self.swap[1]:
+      im1 = im1[:,::-1,:]
+    if roi is None:
+      pass
+    elif isinstance(roi,slice):
+      im1 = im1[:,roi,:]
+      im2 = im2[:,roi,:]
+    elif isinstance(roi,int):
+      if not hasattr(self,"roi1"): self.roi1 = alignment.findRoi(im1[0],roi)
+      if not hasattr(self,"roi2"): self.roi2 = alignment.findRoi(im2[0],roi)
+      im1 = im1[:,self.roi1,:]; im2 = im2[:,self.roi2,:]
+    return im1,im2
+
+  def guiAlign(self,shot=0,save="auto"):
+    im1,im2 = self.getShot(shot)
+    gui = alignment.GuiAlignment(im1[0],im2[0])
+    input("Enter to start")
+    gui.start()
+    if save == "auto":
+      fname = g_folder_init+"/run%04d_gui_align.npy" % self.run
+    else:
+      fname = save
+    self.initAlign = gui.transform
+    gui.save(fname)
+
+
+  def analyzeScan(self,initpars=None,nShotsPerCalib=20,nC=None,doFit=False,fitEveryCalib=False):
+    """ use xhift = None; in this way the fit routine does not try to automatically find the translationx parameter """
+    if initpars is None: initpars= self.initAlign
+    d = self.d
+    if nC is None: nC = d.opal1.nCalib
+    shots = slice(nShotsPerCalib)
+    for i in range(nC):
+      s1,s2 = self.getShot(shots,calib=i)
+      if fitEveryCalib is not False:
+        res = alignment.doShots(s1[:fitEveryCalib],s2[:fitEveryCalib],doFit=True,\
+                initpars=initpars); #
+        idx = np.argmin( [p.fom for p in res] )
+        res = res[idx]
+        initpars = res.final_pars; self.initAlign=res.final_pars
+        print("Calib cycle %d -> %.3f aligned (best FOM: %.2f)" % (i,self.d.scan.scanmotor0_values[i],res.fom))
+      ret = alignment.doShots(s1,s2,initpars=initpars,doFit=doFit)
+      res = alignment.unravel_results(ret)
+      self.results[i] = res
+    return self.results.values()
+
+  def doShot(self,shot=0,calib=None,initpars=None,im1=None,im2=None,doFit=True,show=False,showInit=False,save=False,savePlot="auto"):
+    if initpars is None: initpars= self.initAlign
+    if (im1 is None) or (im2 is None):
+      im1,im2 = self.getShot(shot,calib=calib); im1=im1[0]; im2=im2[0]
+    r = alignment.doShot(im1,im2,initpars,doFit=doFit,show=showInit)
+    im1 = r.im1
+    im2 = r.im2
+    self.initAlign = r.final_pars
+    if show:
+      if savePlot == "auto":
+        savePlot = g_folder_out+"/run%04d_calib%s_shot%04d_fit.png" % (self.run,calib,shot)
+      alignment.plotShot(im1,im2,res=r,save=savePlot)
+    if save: self.saveTransform()
+    return r
+
+  def doShots(self,shots=slice(0,50),calib=None,initpars=None,doFit=False,returnBestTransform=False,unravel=True):
+    if initpars is None: initpars= self.initAlign
+    d = self.d
+    s1,s2 = self.getShot(shots,calib=calib)
+    ret,t = alignment.doShots(s1,s2,initpars=initpars,doFit=doFit,\
+            returnBestTransform=True)
+    if doFit:
+      self.initAlign = t
+    ret_unravel = alignment.unravel_results(ret)
+    # keep it for later !
+    self.results[calib] = ret_unravel
+    if unravel: ret = ret_unravel
+    if returnBestTransform:
+      return ret,t
+    else:
+      return ret
+
+  def save(self,fname="auto",overwrite=False):
+    if fname == "auto":
+      fname = g_folder_out+"/run%04d_analysis.h5" % self.run
+    if os.path.exists(fname) and not overwrite:
+      print("File %s exists, **NOT** saving, use overwrite=True is you want ..."%fname)
+      return
+    if os.path.exists(fname) and overwrite:
+      os.unlink(fname)
+    print("Saving results to %s"%fname)
+    h = h5py.File(fname)
+    h["roi1"] = (self.roi1.start,self.roi1.stop)
+    h["roi2"] = (self.roi2.start,self.roi2.stop)
+    #h["transform"] = self.initAlign
+    for (c,v) in self.results.items():
+      cname = "calib%04d/" % c if isinstance(c,int) else "calib%s/" % c
+      for p,vv in v.items():
+        if p == "parameters":
+          for pname,parray in vv.items():
+            name = cname + p + "/" + pname
+            h[name] = parray
+        else:
+          h[cname + p] = vv
+    h.close()
+
+
+  def saveTransform(self,fname="auto",transform=None):
+    if transform is None: transform = self.initAlign
+    if fname == "auto":
+      fname = g_folder_init+"/run%04d_transform.npy" % self.run
+    print("Saving roi and transformation parameter to %s"%fname)
+    alignment.saveAlignment(fname,self.initAlign,self.roi1,self.roi2)
+
+  def loadTransform(self,fname="auto"):
+    if fname == "auto":
+      fname = g_folder_init+"/run%04d_transform.npy" % self.run
+    elif isinstance(fname,int):
+      fname = g_folder_init+"/run%04d_transform.npy" % fname
+    temp = np.load(fname).item()
+    self.initAlign = temp["transform"]
+    self.roi1 = temp["roi1"]
+    self.roi2 = temp["roi2"]
+    print("init transform and ROIs from %s"%fname)
+
+
+  def clearCache(self):
+    del self.roi1 
+    del self.roi2
+    alignment.clearCache(); # nedded for multiprocessing can leave bad parameters in the cache
+
+  def setDefaultTransform( self ):
+    t = dict( scalex=0.65,rotation=0.0,transx=90, iblur1=4.3,fix_iblur1=False )
+    self.initAlign = t
+    return t
+
+def quick_mec(run,ref=236,divideByRef=False,returnRes=False):
+  """ useful to analyze the runs around 140 (done with the focusing """
+  ref_run = 236
+  h=h5py.File("mecl3616_output/run%04d_analysis.h5" %ref,"r")
+  ref = np.nanmean(h["calibNone"]["ratio"][...],axis=0)
+  r = AnalyzeRun(run,initAlign=ref,swapx=True,swapy=False)
+  res=r.doShots(slice(5),doFit=False)
+  ret = res["ratio"]/ref if divideByRef else res["ratio"]
+  if returnRes:
+    return ret,res
+  else:
+    return ret
+
+def quickAndDirty(run,nShots=300,returnAll=True,doFit=False):
+  """ useful to analyze the runs around 140 (done with the focusing """
+  r = AnalyzeRun(run,swap=True,initAlign=g_folder_init+"/run0144_transform.npy")
+  res=r.doShots(slice(nShots),doFit=doFit)
+  o = alignment.unravel_results(res)
+  ref = np.nanmedian(o["ratio"][:40],0)
+  sam = np.nanmedian(o["ratio"][50:],0)
+  if returnAll:
+    return sam/ref,o["ratio"]/ref
+  else:
+    return sam/ref