Source code for haslib.tof_class

# -*- coding: utf-8 -*-
# tof_class.py
#
# Class definition for an object that should hold all information of an
# inelastic measurement run
#
# Patrick Kraus
# 01.03.2012


from io import StringIO

from matplotlib.pyplot import figure, plot, show, xlabel, ylabel
from numpy import arange, array, exp, hstack, loadtxt, pi, real, where, zeros
from numpy.fft import fft, ifft
from numpy import random
from scipy import interpolate

from .science_basics import fitgauss_cf, kb, mel

# ion() # not sure if correct 


[docs]class TOF(object): # ----------------------------------------------------------------- # ---------------------- INITIALROUTINE --------------------------- # ----------------------------------------------------------------- def __init__(self, appa, G1, G2, DATA, SEQUENCE, specflag=False, accept=False, Mode=1, bounds=[-10.0, 10.0], read_all=True): """ TOF Object - Storing the full information of the TOF run Starting Parameters: appa ... Apparatus Object on which the Measurement was performed G1, G2 ... Reciprocal direction of the measurement run DATA ... filename and adress of the TOF file SEQUENCE ... Entfaltungssequenz (array) specflag ... Boolean defining if the measurement was taken at the specular position. If so, it can be used for nozzle temperature calibration accept ... Toggle if the program should ask for temperaturechanges Mode ... Toggle deconvolution mode 1=fft, everything else uses Huisken-summation """ # Storage of initial parameters self.G = array([G1, G2]) self.appa = appa # Reading the handed file self.file = DATA # Calculating System Time from filename T1, T2 = self.file[self.file.rfind( '_')+1:self.file.rfind('.')].split('.') self.mes_time = float(T1)*1E4+float(T2) # Variable Declaration self.theta_app = None self.delay = None self.slottime = None self.windows = None self.duration = None self.TN = None self.Tsamp = None # Building a container for phonon data self.phoncont = array([0.0, 0.0]) # Variablemap dictionary variable_dict = {'theta /° ': 'theta_app', 'phi /° ': 'theta_app', 'delay /mus ': 'delay', 'slottime /ns ': 'slottime', 'windows': 'windows', 'duration /s': 'duration', 'nozzle temperature in K': 'TN', 'sample temperature in Celsius': 'Tsamp'} HEADER = True WRIT = False hedflag = True self.meas = array([0]) self.toggle = accept # Read the file with open(self.file, encoding='utf-8') as myfile: for line in myfile: arg = line[line.find('[')+1:line.find(']')] if line.find('[') >= 0: if arg in variable_dict: setattr(self, variable_dict[arg], float( line[line.find(']')+1:-1])) if arg == "comment: ": self.TN = float(line[14:18]) else: try: dummy = int(arg[1]) if self.meas.size != 1: if read_all: self.meas = self.meas+loadtxt(StringIO(line)) else: self.meas = loadtxt(StringIO(line)) else: self.meas = loadtxt(StringIO(line)) except ValueError: pass self.specflag = specflag # plot(SEQUENCE) # Nun die Entfaltung des Signals if Mode == 1: self.timesig = real(ifft(fft(SEQUENCE[::-1])*fft(self.meas))) elif Mode == 2: n_deconv = 100.0 # anzahl der beruecksichtigten phasen self.timesig = zeros(SEQUENCE.size) for derun in range(int(n_deconv)): self.timesig = self.timesig + \ real(ifft(fft(SEQUENCE[::-1]) * fft(self.meas)*exp(random.rand()*2*pi*1j))) self.timesig = self.timesig*15/n_deconv else: self.timesig = zeros(self.meas.size) # for ele in range(self.timesig.size): # self.timesig[ele] = 4*sum(SEQUENCE*roll(self.meas,1-ele))/(self.meas.size+1) testsig = hstack([self.meas, self.meas]) for j in range(self.timesig.size): summ = 0 for k in range(self.timesig.size): summ = summ + SEQUENCE[k]*testsig[j+k+1] self.timesig[j] = 4*summ/(self.meas.size+1) # Und die zugehörige Zeit self.time = arange(self.timesig.size)*self.slottime * \ 1E-6 + self.delay*1E-3 # in ms # Nun werden die Parameter der Messapparatur eingelesen self.xcd = self.appa.xcd self.xsd = self.appa.xsd self.xcs = self.appa.xcs # If the specflag is set, the skript automatically calculates the nozzle temperature if self.specflag: self.set_tof_energ() self.appa.specpos = self.theta_app # Now a bool flag is set to read if the signal has already been transformed self.tran = False # Define bounds for the splinefit (should be the same as used in the analysis later) self.lbo = bounds[0] self.ubo = bounds[1] self.rsteps = 1E3 # ----------------------------------------------------------------- # ----------------------------------------------------------------- # ----------------------------------------------------------------- # ----------------------------------------------------------------- # ------------------- MEASUREMENT ROUTINES ------------------------ # -----------------------------------------------------------------
[docs] def set_nozzletemp(self, Tnew): """ Function to set the Nozzle Temperature of the measurement new """ self.TN = Tnew self.appa.set_nozzletemp(Tnew)
[docs] def set_tof_energ(self, std=0.1): """ Function to measure the beam energy used in the TOF measurement. """ if self.specflag: # This is the case when the TOF was made on a specular position # p = fitgauss(self.time,self.timesig-self.timesig.min(),self.time[find(self.timesig==self.timesig.max())][0],self.timesig.max()-self.timesig.min(),0.0001) x0 = self.time[self.timesig.argmax()][0] I = self.timesig.max()-self.timesig.min() y = self.timesig-self.timesig.min() # mask = (self.time > (x0 - std)) & (self.time < (x0 + std)) #mask cutting off in x # mask = y > I/8 #mask cutting off in y # x = self.time[mask] # y = y[mask] # p = fitgauss_cf(x,y,x0,I,std) p = fitgauss_cf(self.time, y, x0, I, std) self.t0 = p[0] v = 1E3*self.xcd/self.t0 Ei = self.appa.m*v**2/2.0 T = 2.0*Ei/(5.0*kb) if self.toggle: self.TN = T else: dummy = input("Change Nozzle Temperature from " + str(self.TN)+" K to "+str(T)+" K ? Y/n >> ") if dummy == 'Y': self.TN = T
# The following proved itself not to be a good idea # self.appa.set_nozzletemp(self.TN)
[docs] def transform_signal(self, mode='new'): """ Function to transform the timesignal to an energy signature """ self.t0 = self.appa.t0 self.dtcd = (self.time-self.t0) self.ttd = self.t0*self.xcs/self.xcd if mode == 'new': self.timesig = self.timesig - \ (self.timesig[-20:-1].mean()+self.timesig[0:19].mean())/2.0 else: self.timesig = self.timesig-self.timesig.min() self.signal_orig = self.timesig * \ ((self.time-self.ttd)*1E-3)**3/(self.xsd**2) self.dE_orig = self.appa.Ei * \ ((1 + self.xcd*self.dtcd/(self.xsd*self.t0))**(-2) - 1)/mel # Perform a splinefit tempcont = array([self.dE_orig, self.signal_orig]) tempcont = tempcont[:, tempcont[0, :].argsort()] self.dE = arange(self.lbo, self.ubo, (self.ubo-self.lbo)/self.rsteps) # ,xb=self.lbo,xe=self.ubo,s=0) tck = interpolate.splrep(tempcont[0, :], tempcont[1, :]) self.signal = interpolate.splev(self.dE, tck) # Now the vectors for constant energy bins # I use a resolution of 0.2 meV binsize = 2 self.dE_eq = arange(self.dE.min(), self.dE.max(), 0.1*binsize) self.signal_eq = zeros(self.dE_eq.size) for index in range(0, self.dE_eq.size): container = array([0]) for bin_i in range(binsize): try: valu = self.dE_eq[index+bin_i].round(1) except IndexError: pass cont = where(self.dE.round(1) == valu)[0] container = hstack([container, cont]) container = container[1:-1] sign = self.signal[container].mean() for bin_i in range(binsize): try: self.signal_eq[index+bin_i] = sign except IndexError: pass # Interpolation Test tck = interpolate.splrep(self.dE_eq, self.signal_eq) self.signal_equ = interpolate.splev(self.dE-0.05, tck) # self.signal_equ = self.signal_eq # Now reset the bool flag self.tran = True
# ----------------------------------------------------------------- # ----------------------------------------------------------------- # ----------------------------------------------------------------- # ----------------------------------------------------------------- # ---------------------- PLOTTING ROUTINES ------------------------ # -----------------------------------------------------------------
[docs] def plot_timesig(self): """ Function to plot the Measured Timesignal """ plot(self.time, self.timesig-self.timesig.min()) xlabel('He travel Time / ms') ylabel('Measured Intensity / a.u.') show()
[docs] def plot_signal(self): """ Function to plot the transformed energy signature """ fig = figure() ax = fig.add_subplot(111) ax.plot(self.dE, self.signal, 'k') ax.set_xlim([self.lbo, self.ubo]) ax.set_ylim([0, 1.1*self.signal.max()]) ax.set_xlabel('Energy Transfer / a.u.') ax.set_ylabel('Intensity / a.u.') show()
# ----------------------------------------------------------------- # ----------------------------------------------------------------- # -----------------------------------------------------------------