import sys
import os
import itertools
import glob
import numpy
import time
#python3 support
PY3 = sys.version_info > (3,)
if PY3:
pass
else:
from itertools import izip as zip
try:
from PyMca import specfilewrapper, EdfFile, SixCircle, specfile
except ImportError:
from PyMca5.PyMca import specfilewrapper, EdfFile, SixCircle, specfile
from .. import backend, errors, util
[docs]class pixels(backend.ProjectionBase):
[docs] def project(self, wavelength, UB, gamma, delta, theta, mu, chi, phi):
y, x = numpy.mgrid[slice(None, gamma.shape[0]), slice(None, delta.shape[0])]
return (y, x)
[docs] def get_axis_labels(self):
return 'y', 'x'
[docs]class HKLProjection(backend.ProjectionBase):
# arrays: gamma, delta
# scalars: theta, mu, chi, phi
[docs] def project(self, wavelength, UB, gamma, delta, theta, mu, chi, phi):
R = SixCircle.getHKL(wavelength, UB, gamma=gamma, delta=delta, theta=theta, mu=mu, chi=chi, phi=phi)
shape = gamma.size, delta.size
H = R[0, :].reshape(shape)
K = R[1, :].reshape(shape)
L = R[2, :].reshape(shape)
return (H, K, L)
[docs] def get_axis_labels(self):
return 'H', 'K', 'L'
[docs]class HKProjection(HKLProjection):
[docs] def project(self, wavelength, UB, gamma, delta, theta, mu, chi, phi):
H, K, L = super(HKProjection, self).project(wavelength, UB, gamma, delta, theta, mu, chi, phi)
return (H, K)
[docs] def get_axis_labels(self):
return 'H', 'K'
[docs]class specularangles(backend.ProjectionBase):
[docs] def project(self, wavelength, UB, gamma, delta, theta, mu, chi, phi):
delta, gamma = numpy.meshgrid(delta, gamma)
mu *= numpy.pi/180
delta *= numpy.pi/180
gamma *= numpy.pi/180
chi *= numpy.pi/180
phi *= numpy.pi/180
theta *= numpy.pi/180
def mat(u, th):
ux, uy, uz = u[0], u[1], u[2]
sint = numpy.sin(th)
cost = numpy.cos(th)
mcost = (1 - numpy.cos(th))
return numpy.matrix([[cost + ux**2 * mcost, ux * uy * mcost - uz * sint, ux * uz * mcost + uy * sint],
[uy * ux * mcost + uz * sint, cost + uy**2 * mcost, uy * uz - ux * sint],
[uz * ux * mcost - uy * sint, uz * uy * mcost + ux * sint, cost + uz**2 * mcost]])
def rot(vx, vy, vz, u, th):
R = mat(u, th)
return R[0, 0] * vx + R[0, 1] * vy + R[0, 2] * vz, R[1, 0] * vx + R[1, 1] * vy + R[1, 2] * vz, R[2, 0] * vx + R[2, 1] * vy + R[2, 2] * vz
#what are the angles of kin and kout in the sample frame?
#angles in the hexapod frame
koutx, kouty, koutz = numpy.sin(- numpy.pi / 2 + gamma) * numpy.cos(delta), numpy.sin(- numpy.pi / 2 + gamma) * numpy.sin(delta), numpy.cos(- numpy.pi / 2 + gamma)
kinx, kiny, kinz = numpy.sin(numpy.pi / 2 - mu), 0, numpy.cos(numpy.pi / 2 - mu)
#now we rotate the frame around hexapod rotation th
xaxis = numpy.array(rot(1, 0, 0, numpy.array([0, 0, 1]), theta))
yaxis = numpy.array(rot(0, 1, 0, numpy.array([0, 0, 1]), theta))
#first we rotate the sample around the xaxis
koutx, kouty, koutz = rot(koutx, kouty, koutz, xaxis, chi)
kinx, kiny, kinz = rot(kinx, kiny, kinz, xaxis, chi)
yaxis = numpy.array(rot(yaxis[0], yaxis[1], yaxis[2], xaxis, chi)) # we also have to rotate the yaxis
#then we rotate the sample around the yaxis
koutx, kouty, koutz = rot(koutx, kouty, koutz, yaxis, phi)
kinx, kiny, kinz = rot(kinx, kiny, kinz, yaxis, phi)
#to calculate the equivalent gamma, delta and mu in the sample frame we rotate the frame around the sample z which is 0,0,1
back = numpy.arctan2(kiny, kinx)
koutx, kouty, koutz = rot(koutx, kouty, koutz, numpy.array([0, 0, 1]), -back)
kinx, kiny, kinz = rot(kinx, kiny, kinz, numpy.array([0, 0, 1]), -back)
mu = numpy.arctan2(kinz, kinx) * numpy.ones_like(delta)
delta = numpy.pi - numpy.arctan2(kouty, koutx)
gamma = numpy.pi - numpy.arctan2(koutz, koutx)
delta[delta > numpy.pi] -= 2 * numpy.pi
gamma[gamma > numpy.pi] -= 2 * numpy.pi
mu *= 1 / numpy.pi * 180
delta *= 1 / numpy.pi * 180
gamma *= 1 / numpy.pi * 180
return (gamma - mu, gamma + mu, delta)
[docs] def get_axis_labels(self):
return 'g-m', 'g+m', 'delta'
[docs]class ThetaLProjection(backend.ProjectionBase):
# arrays: gamma, delta
# scalars: theta, mu, chi, phi
[docs] def project(self, wavelength, UB, gamma, delta, theta, mu, chi, phi):
R = SixCircle.getHKL(wavelength, UB, gamma=gamma, delta=delta, theta=theta, mu=mu, chi=chi, phi=phi)
shape = gamma.size, delta.size
L = R[2, :].reshape(shape)
theta_array = numpy.ones_like(L) * theta
return (theta_array, L)
[docs] def get_axis_labels(self):
return 'Theta', 'L'
[docs]class QProjection(backend.ProjectionBase):
[docs] def project(self, wavelength, UB, gamma, delta, theta, mu, chi, phi):
shape = gamma.size, delta.size
sixc = SixCircle.SixCircle()
sixc.setLambda(wavelength)
sixc.setUB(UB)
R = sixc.getQSurface(gamma=gamma, delta=delta, theta=theta, mu=mu, chi=chi, phi=phi)
qz = R[0, :].reshape(shape)
qy = R[1, :].reshape(shape)
qx = R[2, :].reshape(shape)
return (qz, qy, qx)
[docs] def get_axis_labels(self):
return 'qx', 'qy', 'qz'
[docs]class SphericalQProjection(QProjection):
[docs] def project(self, wavelength, UB, gamma, delta, theta, mu, chi, phi):
qz, qy, qx = super(SphericalQProjection, self).project(wavelength, UB, gamma, delta, theta, mu, chi, phi)
q = numpy.sqrt(qx**2 + qy**2 + qz**2)
theta = numpy.arccos(qz / q)
phi = numpy.arctan2(qy, qx)
return (q, theta, phi)
[docs] def get_axis_labels(self):
return 'Q', 'Theta', 'Phi'
[docs]class CylindricalQProjection(QProjection):
[docs] def project(self, wavelength, UB, gamma, delta, theta, mu, chi, phi):
qz, qy, qx = super(CylindricalQProjection, self).project(wavelength, UB, gamma, delta, theta, mu, chi, phi)
qpar = numpy.sqrt(qx**2 + qy**2)
phi = numpy.arctan2(qy, qx)
return (qpar, qz, phi)
[docs] def get_axis_labels(self):
return 'qpar', 'qz', 'Phi'
[docs]class nrQProjection(QProjection):
[docs] def project(self, wavelength, UB, gamma, delta, theta, mu, chi, phi):
qx, qy, qz = super(nrQProjection, self).project(wavelength, UB, gamma, delta, 0, mu, chi, phi)
return (qx, qy, qz)
[docs] def get_axis_labels(self):
return 'qx', 'qy', 'qz'
[docs]class TwoThetaProjection(SphericalQProjection):
[docs] def project(self, wavelength, UB, gamma, delta, theta, mu, chi, phi):
q, theta, phi = super(TwoThetaProjection, self).project(wavelength, UB, gamma, delta, theta, mu, chi, phi)
return 2 * numpy.arcsin(q * wavelength / (4 * numpy.pi)) / numpy.pi * 180, # note: we need to return a 1-tuple?
[docs] def get_axis_labels(self):
return 'TwoTheta'
[docs]class Qpp(nrQProjection):
[docs] def project(self, wavelength, UB, gamma, delta, theta, mu, chi, phi):
qx, qy, qz = super(Qpp, self).project(wavelength, UB, gamma, delta, theta, mu, chi, phi)
qpar = numpy.sqrt(qx**2 + qy**2)
qpar[numpy.sign(qx) == -1] *= -1
return (qpar, qz)
[docs] def get_axis_labels(self):
return 'Qpar', 'Qz'
[docs]class GammaDeltaTheta(HKLProjection): # just passing on the coordinates, makes it easy to accurately test the theta correction
[docs] def project(self, wavelength, UB, gamma, delta, theta, mu, chi, phi):
delta, gamma = numpy.meshgrid(delta, gamma)
theta = theta * numpy.ones_like(delta)
return (gamma, delta, theta)
[docs] def get_axis_labels(self):
return 'Gamma', 'Delta', 'Theta'
[docs]class GammaDelta(HKLProjection): # just passing on the coordinates, makes it easy to accurately test the theta correction
[docs] def project(self, wavelength, UB, gamma, delta, theta, mu, chi, phi):
delta, gamma = numpy.meshgrid(delta, gamma)
return (gamma, delta)
[docs] def get_axis_labels(self):
return 'Gamma', 'Delta'
[docs]class GammaDeltaMu(HKLProjection): # just passing on the coordinates, makes it easy to accurately test the theta correction
[docs] def project(self, wavelength, UB, gamma, delta, theta, mu, chi, phi):
delta, gamma = numpy.meshgrid(delta, gamma)
mu = mu * numpy.ones_like(delta)
return (gamma, delta, mu)
[docs] def get_axis_labels(self):
return 'Gamma', 'Delta', 'Mu'
[docs]class EH1(ID03Input):
monitor_counter = 'mon'
[docs] def parse_config(self, config):
super(EH1, self).parse_config(config)
self.config.centralpixel = util.parse_tuple(config.pop('centralpixel'), length=2, type=int) # x,y
self.config.hr = config.pop('hr', None) # Optional, hexapod rotations in miliradians. At the entered value the sample is assumed flat, if not entered the sample is assumed flat at the spec values.
self.config.UB = config.pop('ub', None) # Optional, takes specfile matrix by default
if self.config.UB:
self.config.UB = util.parse_tuple(self.config.UB, length=9, type=float)
if self.config.hr:
self.config.hr = util.parse_tuple(self.config.hr, length=2, type=float)
[docs] def process_image(self, scanparams, pointparams, image):
gamma, delta, theta, chi, phi, mu, mon, transm, hrx, hry = pointparams
wavelength, UB = scanparams
weights = numpy.ones_like(image)
if self.config.hr:
zerohrx, zerohry = self.config.hr
chi = (hrx - zerohrx) / numpy.pi * 180. / 1000
phi = (hry - zerohry) / numpy.pi * 180. / 1000
if self.config.background:
data = image / mon
else:
data = image / mon / transm
if mon == 0:
raise errors.BackendError('Monitor is zero, this results in empty output. Scannumber = {0}, pointnumber = {1}. Did you forget to open the shutter?'.format(self.dbg_scanno, self.dbg_pointno))
util.status('{4}| gamma: {0}, delta: {1}, theta: {2}, mu: {3}'.format(gamma, delta, theta, mu, time.ctime(time.time())))
# pixels to angles
pixelsize = numpy.array(self.config.pixelsize)
sdd = self.config.sdd
app = numpy.arctan(pixelsize / sdd) * 180 / numpy.pi
centralpixel = self.config.centralpixel # (column, row) = (delta, gamma)
gamma_range = -app[1] * (numpy.arange(data.shape[1]) - centralpixel[1]) + gamma
delta_range = app[0] * (numpy.arange(data.shape[0]) - centralpixel[0]) + delta
# masking
if self.config.maskmatrix is not None:
if self.config.maskmatrix.shape != data.shape:
raise errors.BackendError('The mask matrix does not have the same shape as the images')
weights *= self.config.maskmatrix
gamma_range = gamma_range[self.config.ymask]
delta_range = delta_range[self.config.xmask]
intensity = self.apply_mask(data, self.config.xmask, self.config.ymask)
weights = self.apply_mask(weights, self.config.xmask, self.config.ymask)
#polarisation correction
delta_grid, gamma_grid = numpy.meshgrid(delta_range, gamma_range)
Pver = 1 - numpy.sin(delta_grid * numpy.pi / 180.)**2 * numpy.cos(gamma_grid * numpy.pi / 180.)**2
intensity /= Pver
return intensity, weights, (wavelength, UB, gamma_range, delta_range, theta, mu, chi, phi)
[docs] def get_point_params(self, scan, first, last):
sl = slice(first, last+1)
GAM, DEL, TH, CHI, PHI, MU, MON, TRANSM, HRX, HRY = list(range(10))
params = numpy.zeros((last - first + 1, 10)) # gamma delta theta chi phi mu mon transm
params[:, CHI] = scan.motorpos('Chi')
params[:, PHI] = scan.motorpos('Phi')
try:
params[:, HRX] = scan.motorpos('hrx')
params[:, HRY] = scan.motorpos('hry')
except:
raise errors.BackendError('The specfile does not accept hrx and hry as a motor label. Have you selected the right hutch? Scannumber = {0}, pointnumber = {1}'.format(self.dbg_scanno, self.dbg_pointno))
if self.is_zap(scan):
if 'th' in scan.alllabels():
th = scan.datacol('th')[sl]
if len(th) > 1:
sign = numpy.sign(th[1] - th[0])
else:
sign = 1
# correction for difference between back and forth in th motor
params[:, TH] = th + sign * self.config.th_offset
else:
params[:, TH] = scan.motorpos('Theta')
params[:, GAM] = scan.motorpos('Gam')
params[:, DEL] = scan.motorpos('Delta')
params[:, MU] = scan.motorpos('Mu')
params[:, MON] = scan.datacol('zap_mon')[sl]
transm = scan.datacol('zap_transm')
transm[-1] = transm[-2] # bug in specfile
params[:, TRANSM] = transm[sl]
else:
if 'hrx' in scan.alllabels():
params[:, HRX] = scan.datacol('hrx')[sl]
if 'hry' in scan.alllabels():
params[:, HRY] = scan.datacol('hry')[sl]
params[:, TH] = scan.datacol('thcnt')[sl]
params[:, GAM] = scan.datacol('gamcnt')[sl]
params[:, DEL] = scan.datacol('delcnt')[sl]
try:
params[:, MON] = scan.datacol(self.monitor_counter)[sl] # differs in EH1/EH2
except:
raise errors.BackendError('The specfile does not accept {2} as a monitor label. Have you selected the right hutch? Scannumber = {0}, pointnumber = {1}'.format(self.dbg_scanno, self.dbg_pointno, self.monitor_counter))
params[:, TRANSM] = scan.datacol('transm')[sl]
params[:, MU] = scan.datacol('mucnt')[sl]
return params
[docs]class EH2(ID03Input):
monitor_counter = 'Monitor'
[docs] def parse_config(self, config):
super(EH2, self).parse_config(config)
self.config.centralpixel = util.parse_tuple(config.pop('centralpixel'), length=2, type=int) # x,y
self.config.UB = config.pop('ub', None) # Optional, takes specfile matrix by default
if self.config.UB:
self.config.UB = util.parse_tuple(self.config.UB, length=9, type=float)
[docs] def process_image(self, scanparams, pointparams, image):
gamma, delta, theta, chi, phi, mu, mon, transm = pointparams
wavelength, UB = scanparams
weights = numpy.ones_like(image)
if self.config.background:
data = image / mon
else:
data = image / mon / transm
if mon == 0:
raise errors.BackendError('Monitor is zero, this results in empty output. Scannumber = {0}, pointnumber = {1}. Did you forget to open the shutter?'.format(self.dbg_scanno, self.dbg_pointno))
util.status('{4}| gamma: {0}, delta: {1}, theta: {2}, mu: {3}'.format(gamma, delta, theta, mu, time.ctime(time.time())))
# area correction
sdd = self.config.sdd / numpy.cos(gamma * numpy.pi / 180)
data *= (self.config.sdd / sdd)**2
# pixels to angles
pixelsize = numpy.array(self.config.pixelsize)
app = numpy.arctan(pixelsize / sdd) * 180 / numpy.pi
centralpixel = self.config.centralpixel # (row, column) = (gamma, delta)
gamma_range = - 1 * app[0] * (numpy.arange(data.shape[0]) - centralpixel[0]) + gamma
delta_range = app[1] * (numpy.arange(data.shape[1]) - centralpixel[1]) + delta
# masking
if self.config.maskmatrix is not None:
if self.config.maskmatrix.shape != data.shape:
raise errors.BackendError('The mask matrix does not have the same shape as the images')
weights *= self.config.maskmatrix
gamma_range = gamma_range[self.config.xmask]
delta_range = delta_range[self.config.ymask]
intensity = self.apply_mask(data, self.config.xmask, self.config.ymask)
weights = self.apply_mask(weights, self.config.xmask, self.config.ymask)
intensity = numpy.fliplr(intensity)
intensity = numpy.rot90(intensity)
weights = numpy.fliplr(weights) # TODO: should be done more efficiently. Will prob change with new HKL calculations
weights = numpy.rot90(weights)
#polarisation correction
delta_grid, gamma_grid = numpy.meshgrid(delta_range, gamma_range)
Phor = 1 - (numpy.sin(mu * numpy.pi / 180.) * numpy.sin(delta_grid * numpy.pi / 180.) * numpy.cos(gamma_grid * numpy.pi / 180.) + numpy.cos(mu * numpy.pi / 180.) * numpy.sin(gamma_grid * numpy.pi / 180.))**2
intensity /= Phor
return intensity, weights, (wavelength, UB, gamma_range, delta_range, theta, mu, chi, phi)
[docs] def get_point_params(self, scan, first, last):
sl = slice(first, last+1)
GAM, DEL, TH, CHI, PHI, MU, MON, TRANSM = list(range(8))
params = numpy.zeros((last - first + 1, 8)) # gamma delta theta chi phi mu mon transm
params[:, CHI] = scan.motorpos('Chi')
params[:, PHI] = scan.motorpos('Phi')
if self.is_zap(scan):
if 'th' in scan.alllabels():
th = scan.datacol('th')[sl]
if len(th) > 1:
sign = numpy.sign(th[1] - th[0])
else:
sign = 1
# correction for difference between back and forth in th motor
params[:, TH] = th + sign * self.config.th_offset
else:
params[:, TH] = scan.motorpos('Theta')
params[:, GAM] = scan.motorpos('Gamma')
params[:, DEL] = scan.motorpos('Delta')
params[:, MU] = scan.motorpos('Mu')
params[:, MON] = scan.datacol('zap_mon')[sl]
transm = scan.datacol('zap_transm')
transm[-1] = transm[-2] # bug in specfile
params[:, TRANSM] = transm[sl]
else:
params[:, TH] = scan.datacol('thcnt')[sl]
params[:, GAM] = scan.datacol('gamcnt')[sl]
params[:, DEL] = scan.datacol('delcnt')[sl]
try:
params[:, MON] = scan.datacol(self.monitor_counter)[sl] # differs in EH1/EH2
except:
raise errors.BackendError('The specfile does not accept {2} as a monitor label. Have you selected the right hutch? Scannumber = {0}, pointnumber = {1}'.format(self.dbg_scanno, self.dbg_pointno, self.monitor_counter))
params[:, TRANSM] = scan.datacol('transm')[sl]
params[:, MU] = scan.datacol('mucnt')[sl]
return params
[docs]class GisaxsDetector(ID03Input):
monitor_counter = 'mon'
[docs] def process_image(self, scanparams, pointparams, image):
ccdy, ccdz, theta, chi, phi, mu, mon, transm = pointparams
weights = numpy.ones_like(image)
wavelength, UB = scanparams
if self.config.background:
data = image / mon
else:
data = image / mon / transm
if mon == 0:
raise errors.BackendError('Monitor is zero, this results in empty output. Scannumber = {0}, pointnumber = {1}. Did you forget to open the shutter?'.format(self.dbg_scanno, self.dbg_pointno))
util.status('{4}| ccdy: {0}, ccdz: {1}, theta: {2}, mu: {3}'.format(ccdy, ccdz, theta, mu, time.ctime(time.time())))
# pixels to angles
pixelsize = numpy.array(self.config.pixelsize)
sdd = self.config.sdd
directbeam = (self.config.directbeam[0] - (ccdy - self.config.directbeam_coords[0]) * pixelsize[0], self.config.directbeam[1] - (ccdz - self.config.directbeam_coords[1]) * pixelsize[1])
gamma_distance = - pixelsize[1] * (numpy.arange(data.shape[1]) - directbeam[1])
delta_distance = - pixelsize[0] * (numpy.arange(data.shape[0]) - directbeam[0])
gamma_range = numpy.arctan2(gamma_distance, sdd) / numpy.pi * 180 - mu
delta_range = numpy.arctan2(delta_distance, sdd) / numpy.pi * 180
#sample pixel distance
spd = numpy.sqrt(gamma_distance**2 + delta_distance**2 + sdd**2)
data *= spd**2 / sdd
# masking
if self.config.maskmatrix is not None:
if self.config.maskmatrix.shape != data.shape:
raise errors.BackendError('The mask matrix does not have the same shape as the images')
weights *= self.config.maskmatrix
gamma_range = gamma_range[self.config.ymask]
delta_range = delta_range[self.config.xmask]
intensity = self.apply_mask(data, self.config.xmask, self.config.ymask)
weights = self.apply_mask(weights, self.config.xmask, self.config.ymask)
return intensity, weights, (wavelength, UB, gamma_range, delta_range, theta, mu, chi, phi)
[docs] def parse_config(self, config):
super(GisaxsDetector, self).parse_config(config)
self.config.directbeam = util.parse_tuple(config.pop('directbeam'), length=2, type=int)
self.config.directbeam_coords = util.parse_tuple(config.pop('directbeam_coords'), length=2, type=float) # Coordinates of ccdy and ccdz at the direct beam position
[docs] def get_point_params(self, scan, first, last):
sl = slice(first, last+1)
CCDY, CCDZ, TH, CHI, PHI, MU, MON, TRANSM = list(range(8))
params = numpy.zeros((last - first + 1, 8)) # gamma delta theta chi phi mu mon transm
params[:, CHI] = scan.motorpos('Chi')
params[:, PHI] = scan.motorpos('Phi')
params[:, CCDY] = scan.motorpos('ccdy')
params[:, CCDZ] = scan.motorpos('ccdz')
params[:, TH] = scan.datacol('thcnt')[sl]
try:
params[:, MON] = scan.datacol(self.monitor_counter)[sl] # differs in EH1/EH2
except:
raise errors.BackendError('The specfile does not accept {2} as a monitor label. Have you selected the right hutch? Scannumber = {0}, pointnumber = {1}'.format(self.dbg_scanno, self.dbg_pointno, self.monitor_counter))
params[:, TRANSM] = scan.datacol('transm')[sl]
params[:, MU] = scan.datacol('mucnt')[sl]
return params
[docs] def find_edfs(self, pattern, scanno):
files = glob.glob(pattern)
ret = {}
for file in files:
try:
filename = os.path.basename(file).split('.')[0]
scan, point = filename.split('_')[-2:]
scan, point = int(scan), int(point)
if scan == scanno and point not in list(ret.keys()):
ret[point] = file
except ValueError:
continue
return ret
[docs]def load_matrix(filename):
if filename == None:
return None
if os.path.exists(filename):
ext = os.path.splitext(filename)[-1]
if ext == '.txt':
return numpy.array(numpy.loadtxt(filename), dtype=numpy.bool)
elif ext == '.npy':
return numpy.array(numpy.load(filename), dtype=numpy.bool)
elif ext == '.edf':
return numpy.array(EdfFile.EdfFile(filename).getData(0), dtype=numpy.bool)
else:
raise ValueError('unknown extension {0}, unable to load matrix!\n'.format(ext))
else:
raise IOError('filename: {0} does not exist. Can not load matrix'.format(filename))