Source code for haslib.apparatus_class

# -*- 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()
# ----------------------------------------------------------------- # ----------------------------------------------------------------- # -----------------------------------------------------------------