# -*- coding: utf-8 -*-
# apparatus_class.py
#
# Class definition for an object which simulates the measurement
# apparatus and its parameters
#
# Patrick Kraus
# 27.02.2012
# import the used libraries
from __future__ import division, print_function
import os
from haslib.science_basics import (mel, m_He4, kb, hb, meV, find_inters,
myunique, gauss, fitgauss, fitgigagauss,
gigagauss)
from numpy import (pi, sqrt, array, sin, arange, vstack,
arcsin, nan, cos, linspace, mean, zeros, hstack)
from numpy.linalg import norm
from scipy.optimize import fsolve
from matplotlib.pyplot import (figure, show, ion, rc, title, cm, close,
subplot2grid)
# from numpy import nanmin, nanmax
from PIL import Image
ion()
[docs]class Apparatus(object):
# -----------------------------------------------------------------
# ---------------------- INITIALROUTINE ---------------------------
# -----------------------------------------------------------------
def __init__(self, sample, nozzletemp=75.0, TSD=91.5, xcd=1.643,
xcs=0.529, xsd=1.114, m=m_He4):
"""
apparatus object
Args:
sample: crystal object defining the sample
nozzletemp: Nozzle Temperature
TSD: Fixed apparative Angle (default 91.5 degrees)
xcd: Distance in between the Chopper and the detector
(default 1.643 m)
xcs: Distance in between the Chopper and the source
(default 0.529 m) NOT used in class
xcd: Distance in between the source and the detector
(default 1.114 m) NOT used in class
m: Atom mass of the beam
"""
# Defining the specific variables
self.sample = sample
self.TSD = TSD
self.specpos = 0.0
self.m = m # kg (Helium Mass)
self.TN = nozzletemp # K Nozzle Temperature
self.Ei = 5 * kb * self.TN / 2.0 # Helium Atom Energy in J
self.ki = sqrt(2 * self.m * self.Ei) / hb # Helium momentum m**-1
self.v = sqrt(2 * self.Ei / self.m) # Helium velocity
self.t0 = 1E3 * xcd / self.v
# Apparative constants for TOF mesurements
self.xcd = xcd
self.xcs = xcs
self.xsd = xsd
# Constants for TOF plot
self.tof_lbo = -10.0
self.tof_ubo = 10.0
# Container for phonon data
self.phoncont = array([0.0, 1.0])
self.gausscont = array([0.0, 0.0])
# -----------------------------------------------------------------
# -----------------------------------------------------------------
# -----------------------------------------------------------------
# -----------------------------------------------------------------
# ---------------- ACTUALIZATION - ROUTINES -----------------------
# -----------------------------------------------------------------
[docs] def set_nozzletemp(self, temp):
"""Function to set a new Nozzle Temperature"""
self.TN = temp # K Nozzle Temperature
self.Ei = 5 * kb * self.TN / 2.0 # Helium Atom Energy in J
self.ki = sqrt(2 * self.m * self.Ei) / hb # Helium momentum m**-1
self.v = sqrt(2 * self.Ei / self.m) # Helium velocity
self.t0 = 1E3 * self.xcd / self.v
[docs] def set_specpos(self, newpos):
"""Function to set the apparative specular position"""
self.specpos = newpos
# -----------------------------------------------------------------
# -----------------------------------------------------------------
# -----------------------------------------------------------------
# -----------------------------------------------------------------
# ----------------- PREDICTION - ROUTINES -------------------------
# -----------------------------------------------------------------
[docs] def predict_boundstates(self, angles, direction):
"""
Calculate boundstate energies
Args:
angles: List of angles where SARs are found in degree
direction: Tuple of two indexes in the direction of
the measurement (e.g. (1, 1))
Yields:
scalar: bound state energy in meV
"""
g1, g2 = direction
self.sample.rotate_crystal(g1, g2)
for theta in angles:
theta = array(theta)
for g in self.sample.Gvecs:
curr_g = self.sample.givevec(g[0], g[1])
bse = ((self.ki * sin(theta * pi / 180.0) + curr_g[0]) ** 2 +
(curr_g[1]) ** 2 - self.ki ** 2) * \
(hb ** 2 / (2 * self.m))
if bse < 0:
bse = nan
yield bse / mel
[docs] def sim_elastic_measurement(self, G1, G2, resflag=False, plotit=True,
fullout=False):
"""
Simulate an elastic measurement in the <G1, G2> direction.
Args:
G1, G2: Direction in which the measurement should
be performed
resflag: Boolean Flag to turn the calculation of resonance
positions on and off
plotit: Boolean Flag to turn the plot off
fullout: Boolean Flag to turn the output of expected angles on
"""
# Turning the crystal in the desired direction
self.sample.rotate_crystal(G1, G2)
# Now the calculation of the positions of the elastic scattering peaks
elastics = array([0, 0, 0, 0])
ASD = self.TSD * pi / 180.0
Gnow = self.sample.givevec(G1, G2)
for iterator in arange(-15, 16, 1):
# removed minus in front of iterator. I think it's wrong - Michael
# angle = 180 * (ASD / 2.0 - arcsin(iterator * (Gnow[0]) /
# (2 * self.ki * cos(ASD / 2.0)))) / pi
# changed above equation to the below one - I think this
# should be norm(G) instead of G(1) - Toni
angle = 180*(ASD/2.0-arcsin(iterator*(norm(Gnow)) /
(2 * self.ki * cos(ASD / 2.0)))) / pi
# print(angle)
if (angle > -16) & (angle < 90):
elastics = vstack([elastics,
[angle, -1, iterator * G1, iterator * G2]])
elastics = elastics[1::]
# Now for the prediction of the resonance positions
rescontainer = array([0, 0, 0, 0])
for reind in range(self.sample.interactE.size):
for gind in range(self.sample.Gvecs.shape[0]):
# I will need abbrivations
rese = self.sample.interactE[reind]
Ga = self.sample.Gvecs[gind, 0]
Gb = self.sample.Gvecs[gind, 1]
# Now I should find this vectors components in the current
# orientation
currG = self.sample.givevec(Ga, Gb)
# Now model the angle-energy resonance function minus the BSE
def res_ang_energ(anglein):
return (((self.ki * sin(anglein * pi / 180.0) +
currG[0]) ** 2 +
(currG[1]) ** 2 - self.ki ** 2) *
(hb ** 2 / (2 * self.m)) - abs(rese))
# Now find the zero-position if there is any
y, dict, du, my = fsolve(res_ang_energ, 0.1,
full_output=1, maxfev=1000)
try:
y = y[0]
except IndexError:
y = y
if (du == 1) & (y > 0) & (y < 90):
rescontainer = vstack([rescontainer, [y, reind, Ga, Gb]])
rescontainer = rescontainer[1::]
# Here also kinematical focussing effects could be calculated
# Everything immediately necessary is calculated right now -
# now for the plotting
if plotit:
simel = figure()
elax = simel.add_subplot(111)
recoldict = {0: 'm', 1: 'b', 2: 'r', 3: 'g', 4: 'y',
5: 'c', 6: '# FF6600', 7: '# ff69b4', -1: 'k'}
for pos in range(elastics.shape[0]):
elax.plot([elastics[pos, 0], elastics[pos, 0]], [0, 10], 'k')
if resflag:
for pos in range(rescontainer.shape[0]):
elax.plot([rescontainer[pos, 0], rescontainer[pos, 0]],
[0, 10], recoldict[int(rescontainer[pos, 1])])
elax.set_xlabel('Incident Angle / Degrees')
elax.set_ylabel('Measured Intensity / cps')
show()
if fullout:
return vstack([elastics, rescontainer])
[docs] def sim_tof_measurement(self, G1, G2, angin, resflag=True, plotit=False,
fullout=True):
"""
Simulate inelastic TOF measurement in the <G1, G2> direction
Args:
resflag: Boolean flag to turn the calculation of resonance
positions on and off. Without this, only the
expected position of the diffuse elastic peak
is calculated (which is always at zero)
G1, G2: Reciprocal orientation of the crystal
angin: Angle of the incoming beam
plotit: Boolean flag to turn plotting on and off
fullout: Boolean flag to turn
"""
# Turning the crystal in the desired direction
self.sample.rotate_crystal(G1, G2)
# Defining output container (also ED peak)
container = array([0.0, -1, 0, 0])
# Calculating resonance positions, if the flag is set
if resflag:
# Loop over the interaction energies:
for reind in range(self.sample.interactE.size):
# Loop over relevant reciprocal lattice vectors
for gind in range(self.sample.Gvecs.shape[0]):
# I will need abbrivations
rese = self.sample.interactE[reind]
Ga = self.sample.Gvecs[gind, 0]
Gb = self.sample.Gvecs[gind, 1]
# Now I should find this vectors components in the current
# orientation
currG = self.sample.givevec(Ga, Gb)
# Now define the scan curve and the resonance curve
def scanc(Q):
return (hb**2*((self.ki*sin(angin*pi/180.0)+Q)**2 /
(sin((self.TSD-angin)*pi/180.0)**2) -
self.ki**2)/(2.0*self.m*mel))
# now the resonance curve
def curveER(Q):
return ((-self.Ei+hb**2*(self.ki*sin(angin*pi/180.0) +
currG[0]+Q)**2/(2.0*self.m)+hb**2*currG[1]**2 /
(2.0*self.m)-rese)/mel)
# Find the intersection (if there is any)
y, du = find_inters(curveER, scanc, 0.1 * currG[0])
# print y
try:
y = y[0]
except IndexError:
y = y
# y = y[0]
if ((du == 1) & (scanc(y) > self.tof_lbo) &
(scanc(y) < self.tof_ubo)):
container = vstack([container,
[scanc(y), reind, Ga, Gb]])
# Remove equal or similarily equal entries
container = myunique(container)
# Now the relevant positions should be calculated, now for the plotting
# and the returns
if plotit:
simtof = figure()
tofax = simtof.add_subplot(111)
recoldict = {0: 'm', 1: 'b', 2: 'r', 3: 'g', 4: 'y', 5: 'c',
6: '# FF6600', 7: '# ff69b4', -1: 'k'}
if container.ndim > 1:
for pos in range(container.shape[0]):
tofax.plot([container[pos, 0], container[pos, 0]],
[0, 10], recoldict[int(container[pos, 1])])
else:
tofax.plot([container[0], container[0]],
[0, 10], recoldict[int(container[1])])
tofax.set_xlabel('Energy Transfer / meV')
tofax.set_ylabel('Intensity / a.u.')
tofax.set_xlim([self.tof_lbo, self.tof_ubo])
show()
if fullout:
return container
# -----------------------------------------------------------------
# -----------------------------------------------------------------
# -----------------------------------------------------------------
# -----------------------------------------------------------------
# ----------------- ANALYSATION - ROUTINES -------------------------
# -----------------------------------------------------------------
[docs] def analyze_elastic(self, El_obj, zoom=10.0, resplflag=True,
useaxes=None, elastic_only=False, fontsize=10):
"""
Function to plot the outcome of an elastic measurement with the
calculated resonance positions
"""
rc('text', usetex=True)
rc('font', family='serif')
# First change the nozzle temperature andf the orientation to
# the values given in El_obj
self.set_nozzletemp(El_obj.TN)
self.sample.rotate_crystal(El_obj.G[0], El_obj.G[1])
# Now calculate the suspected elastics and resonance
# positions for this measurement
POINTS = self.sim_elastic_measurement(El_obj.G[0], El_obj.G[1],
True, False, True)
# Now the Plot routine
colors = cm.rainbow(linspace(0, 1, len(self.sample.interactE) + 1))
if useaxes is None:
mesel = figure()
mesax = mesel.add_subplot(111)
else:
mesax = useaxes
# make title
if(El_obj.G[0] == 1 and El_obj.G[1] == 1):
str_dir = "$\overline{\Gamma \mathrm{K}}$"
else:
str_dir = "$\overline{\Gamma \mathrm{M}}$"
tex_filename = os.path.basename(El_obj.file).replace("_", "\_")
title(("${name}$ @ ${temp:0.1f}^\circ$C," +
"{direction}: T$_N$={nozzle:0.2f}K~({energy:0.2f}meV)," +
"t={dwell:0.2f} s\n" +
" file: {filename}").format(name=self.sample.name,
temp=El_obj.TS,
nozzle=El_obj.TN,
dwell=mean(El_obj.time),
direction=str_dir,
energy=self.Ei/meV,
filename=tex_filename),
fontsize=1.2*fontsize)
mesax.tick_params(axis='both', which='major', labelsize=fontsize)
# Plot error band
cps_std = sqrt(El_obj.counts)/El_obj.time
mesax.fill_between(El_obj.theta + self.TSD / 2.0, El_obj.cps + cps_std,
El_obj.cps - cps_std, facecolor='lightgray')
# Plot the signal
mesax.plot(El_obj.theta + self.TSD / 2.0, El_obj.cps, 'k.-',
linewidth=1, markersize=5)
# title(os.path.basename(El_obj.file), fontsize=fontsize*1.2)
maxtheta = El_obj.theta.max() + self.TSD / 2.0
mintheta = El_obj.theta.min() + self.TSD / 2.0
# Plot the lines
if resplflag:
for pos in range(POINTS.shape[0]):
if POINTS[pos, 0] > maxtheta:
continue
if POINTS[pos, 0] < mintheta:
continue
if elastic_only and POINTS[pos, 1] != -1:
continue
c = colors[int(POINTS[pos, 1])]
mesax.axvline(x=POINTS[pos, 0], color=c, linewidth=1)
n_bse = int(POINTS[pos, 1] - 1)
g1 = int(POINTS[pos, 2])
g2 = int(POINTS[pos, 3])
if POINTS[pos, 1] == -1:
labtext = 'EL\_({:d}, {:d})'.format(g1, g2)
else:
labtext = '({:d}, {:d})\_{:d}'.format(g1, g2, n_bse)
mesax.text(POINTS[pos, 0]+0.2, El_obj.cps.max()*1.05/zoom,
labtext, rotation='vertical', color=c)
# Now for the zoom
mesax.set_xlim([self.TSD / 2.0 + El_obj.theta.min(),
self.TSD / 2.0 + El_obj.theta.max()])
mesax.set_ylim([El_obj.cps.min(), El_obj.cps.max() * 1.1 / zoom])
mesax.set_xlabel('Incident Angle / $^\circ$', fontsize=fontsize)
mesax.set_ylabel('Measured Intensity / a.u.', fontsize=fontsize)
return mesax
[docs] def analyze_tof(self, Tof_obj, ret_fig=False, upper=0, textpos=0,
useaxes=None, mode='new'):
"""
Function to plot the TOF measurement
with the calculated resonance positions
"""
# First check if the TOF object is sufficiently transformed
if Tof_obj.tran:
# Now predict the resonance positions
cont = self.sim_tof_measurement(Tof_obj.G[0], Tof_obj.G[1],
self.TSD / 2.0 -
(Tof_obj.theta_app -
self.specpos))
# Now the plot routine
if useaxes is None:
toffig = figure()
tofax = toffig.add_subplot(111)
else:
tofax = useaxes
# Plot the transformed signal
if mode == 'new':
tofax.plot(Tof_obj.dE, Tof_obj.signal, color='0.75', picker=5)
tofax.plot(Tof_obj.dE, Tof_obj.signal_equ, 'k')
else:
tofax.plot(Tof_obj.dE, Tof_obj.signal, 'k')
self.lbo = 0
if upper == 0:
self.ubo = 1.1 * Tof_obj.signal[0:750].max()
else:
self.ubo = upper
if textpos == 0:
rang = 0.85 * self.ubo - self.lbo
else:
rang = textpos * self.ubo - self.lbo
# Now plot the ED and the resonances
recoldict = {0: 'm', 1: 'b', 2: 'r', 3: 'g', 4: 'y', 5: 'c',
6: '# FF6600', 7: '# ff69b4', -1: 'k'}
if cont.ndim > 1:
for pos in range(cont.shape[0]):
tofax.plot([cont[pos, 0], cont[pos, 0]],
[self.lbo, self.ubo],
recoldict[int(cont[pos, 1])])
if cont[pos, 1] == -1:
labtext = 'EL_(' + \
str(int(cont[pos, 2])) + \
', ' + \
str(int(cont[pos, 3])) + \
')'
elif cont[pos, 1] == 0:
labtext = '(' + str(int(cont[pos, 2])) + \
', ' + str(int(cont[pos, 3])) + ')_T'
else:
labtext = '(' + str(int(cont[pos, 2])) + \
', ' + str(int(cont[pos, 3])) + ')_' + \
str(int(cont[pos, 1] - 1))
tofax.text(cont[pos, 0] + 0.1, rang + self.lbo,
labtext, rotation='vertical',
color=recoldict[int(cont[pos, 1])])
else:
tofax.plot([cont[0], cont[0]], [self.lbo, self.ubo],
recoldict[int(cont[1])])
tofax.set_xlim([self.tof_lbo, self.tof_ubo])
tofax.set_ylim([self.lbo, self.ubo])
tofax.set_xlabel('Energy Transfer / meV')
tofax.set_ylabel('Intensity / a.u.')
else:
print("Object not sufficiently transformed." +
" Run transform_signal()")
if ret_fig:
return toffig
[docs] def prettyprint_analyze_tof(self, TOF_obj, path_to_structuregraph=None,
ylim=30, samplename='Bi$_2$Se$_3$(111)'):
"""
Function to plot the TOF measurements according
to Prof. Benedek's wishes
"""
fig = figure()
ax_tex = subplot2grid((4, 4), (0, 0), rowspan=1, colspan=3)
ax_pl = subplot2grid((4, 4), (1, 0), rowspan=3, colspan=4)
self.analyze_tof(TOF_obj, upper=ylim, useaxes=ax_pl, mode='new')
if path_to_structuregraph is not None:
ax_dir = subplot2grid((4, 4), (0, 3), rowspan=1, colspan=1)
Img = Image.open(path_to_structuregraph)
ax_dir.imshow(Img)
ax_dir.get_xaxis().set_visible(False)
ax_dir.get_yaxis().set_visible(False)
ax_dir.axis('off')
Filname = TOF_obj.file[-1 * TOF_obj.file[::-1].find('/'):-1]
Filname += TOF_obj.file[-1]
Ang = (self.TSD / 2.0 - (TOF_obj.theta_app - self.specpos))
Ei = 5 * kb * self.TN / 2.0
ax_tex.text(0, 0,
samplename+' \n' + Filname + '\n' + '\n' +
r'$E_i$ = ' + str(Ei / mel)[0:5] + ' meV \n' +
r'$\theta_i$ = ' + str(Ang) + r'$^\circ$' + '\n' +
'Sample Temperature = ' +
str(273.15 + TOF_obj.Tsamp)[0:5] +
' K\n' + r'$\theta_{SD}$ = ' +
str(self.TSD) + r'$^\circ$')
ax_tex.get_xaxis().set_visible(False)
ax_tex.get_yaxis().set_visible(False)
ax_tex.axis('off')
return fig
# -----------------------------------------------------------------
# -----------------------------------------------------------------
# -----------------------------------------------------------------
# -----------------------------------------------------------------
# ----------------- MEASUREMENT - ROUTINES ------------------------
# -----------------------------------------------------------------
[docs] def phon_pick_event(self, event):
thisline = event.artist
xdata, ydata = thisline.get_data()
ind = event.ind[0]
self.phoncont = vstack([self.phoncont, array([xdata[ind],
ydata[ind]])])
ax = thisline.get_axes()
ax.plot(xdata[ind], ydata[ind], 'bo')
ax.set_xlim([self.tof_lbo, self.tof_ubo])
ax.set_ylim([self.lbo, self.ubo])
show()
[docs] def gauss_pick_event(self, event):
thisline = event.artist
xdata, ydata = thisline.get_data()
EN = xdata[ydata.argmax()][0]
# Now fit a gaussian to this gaussian
p_sl = fitgauss(xdata, ydata, EN, ydata.max())
self.gausscont = vstack([self.gausscont, array([EN, p_sl[2] ** 2])])
ind = event.ind[0]
ax = thisline.get_axes()
ax.plot(xdata[ind], ydata[ind], 'ro')
ax.set_xlim([self.tof_lbo, self.tof_ubo])
ax.set_ylim([self.lbo, self.ubo])
show()
[docs] def phonlog(self, Tofobj):
"""
Measurement routine to fit and log the positions of phonon
peaks in a TOF energy signature. The given Object must be a
transformed TOF measurement.
"""
# First we check if the conditions for a successfull measurement
# are fulfilled
if Tofobj.tran is False:
print('Error: Object is not sufficiently transformed')
# Write an explanation
print('''
Welcome to phonlog!
Click on the positions where the script should fit a gauss-function
(the ED peak is fitted automatically)
When you are done, write 'done' and press Enter''')
gaussflag = True
while gaussflag:
# Now we perform a standard TOF analysis - and keep the figure
# handle
ion()
toffig = self.analyze_tof(Tofobj, True)
# Initiate the input loop and log the positions
self.phoncont = array([0.0, 1.0])
self.gausscont = array([0.0, 0.0])
cid = toffig.canvas.mpl_connect('pick_event', self.phon_pick_event)
Flaggy = True
while Flaggy:
toffig.canvas.draw()
dummy = input('>>> ')
if dummy == 'done':
Flaggy = False
toffig.canvas.mpl_disconnect(cid)
# Now we need fit the plotted area with the gigagauss function
# First, we prepare the parameters vector
p0 = array([0.0, 0.0, 2.0])
self.p = zeros(p0.shape)
# Now for every chosen point one gauss parameters list
for index in arange(self.phoncont.shape[0]) - 1:
p0 = hstack([p0, array([self.phoncont[index + 1, 0],
self.phoncont[index + 1, 1],
0.1])])
self.p, self.cov = fitgigagauss(Tofobj.dE, Tofobj.signal, p0)
close()
toffig = self.analyze_tof(Tofobj, True)
# Now plot the fitted gaussians
ax = toffig.get_axes()[0]
xx = Tofobj.dE
ax.plot(Tofobj.dE, gauss(array([0, self.p[1], 0.1]), xx), 'k')
ax.plot(Tofobj.dE, gigagauss(self.p, xx), 'r')
for ind in arange((self.p.size - 3) / 3) + 1:
ing = ind * 3
ax.plot(xx, gauss(array([self.p[ing],
abs(self.p[ing + 1]),
self.p[ing + 2]]), xx),
'g', picker=5)
ax.set_xlim([self.tof_lbo, self.tof_ubo])
ax.set_ylim([self.lbo, self.ubo])
print('Enter Y for OK, N to start again')
dummy = input('>>> ')
if dummy == 'Y':
gaussflag = False
else:
close()
# Now the gausspeaks are fitted. The user should now be able to
# select single gausspeaks to save the corresponding positions
# in the TOF file
# Now we initiate the mapping event
print("""
Please mark the gaussians that should be regarded, then write 'done'
""")
cid = toffig.canvas.mpl_connect('pick_event', self.gauss_pick_event)
Gauflag = True
while Gauflag:
dummy = input('>>> ')
if dummy == 'done':
Gauflag = False
toffig.canvas.mpl_disconnect(cid)
print("are these points ok? (Y to accept, N to terminate)")
dummy = input('>>> ')
if dummy == 'Y':
Tofobj.phoncont = self.gausscont[1:, :]
print('logged')
else:
print("Terminating")
self.phoncont = array([0.0, 0.0])
close()
# -----------------------------------------------------------------
# -----------------------------------------------------------------
# -----------------------------------------------------------------