Source code for haslib.corrugation

from __future__ import division
import sys
import os
from haslib.crystal_class import Crystal
from haslib.science_basics import A, meV, h2m, kb, hb
import numpy as np
from haslib.icc_2d_2 import iCC_2D_2

# constants
GM = (1, 0)
GK = (1, 1)


# data classes
[docs]class Peak: def __init__(self, angle, g1, g2, area, err_area): self.angle = angle self.g1 = g1 self.g2 = g2 self.area = area self.err_area = err_area
[docs]class Spectrum: def __init__(self, direction, energy, T, a=None): self.direction = direction self.energy = energy self.T = T self.peaks = [] self.a = a
[docs] def addPeak(self, angle, g1, g2, area, err_area): self.peaks.append(Peak(angle, g1, g2, area, err_area))
[docs]def get_alpha(D, E, theta, slope, m, TSD): """ calculate Td**2 * M from Deybe Waller measurement for a given D """ kib = np.sqrt(E/h2m(m)) kzi = kib * np.sqrt(np.cos(theta*np.pi/180.0)**2 + D/E) # D/E = beeby kzf = kib * np.sqrt(np.cos((TSD-theta)*np.pi/180.0)**2 + D/E) alpha = -(3.0 * hb**2 * (kzi + kzf)**2) / (kb * slope) # = 2W return alpha
[docs]def signal_per_spectrum(spectrum, D, xi, h, a, z, gclose, alpha, m, TSD, st_p_wl): """Calculate signal for each peak in a spectrum Args: spectrum: Spectrum class D: Depth of the Potential xi: Potential parameter h: corrugation a: lattice constant z: (min, max) gclose: closed channels cc alpha: Mass_sample * debye temp ** 2 m: mass of the beam TSD: angle between source and detector st_p_wl: steps per wavelength Returns: array: array each row containing the g1, g2, corrected signal and peak.area """ # check if the spektrum has his own lattice constant if spectrum.a: a = a else: a = spectrum.a # create Crystal class and extract g1, g2 sample = Crystal(a, a, 120) sample.rotate_crystal(*spectrum.direction) g1 = sample.b1 g2 = sample.b2 # alias values from spectrum class Ei = spectrum.energy*meV TS = spectrum.T # use all open channels gopen = -1 # allocate array for g1, g2, corrected signal, peak area signals = np.zeros([len(spectrum.peaks), 4]) kib = np.sqrt(Ei/h2m(m)) for index, peak in enumerate(spectrum.peaks): theta = peak.angle # simulate peak sig = iCC_2D_2(h, 0, a, D, xi, Ei, theta, g1, g2, gopen, gclose, m=m, floq=0, phon=[0, 0, 0], sta_z=z[0], end_z=z[1], st_p_wl=st_p_wl)[:, [0, 1, 6, 5]] # get peak specified in spectrum sig = sig[(sig[:, 0] == int(peak.g1)) & (sig[:, 1] == int(peak.g2)), :] # calculate Debye Waller correction kzi = kib * np.sqrt(np.cos(theta*np.pi/180.0)**2 + D/Ei) kzf = kib * np.sqrt(np.cos((TSD-theta)*np.pi/180.0)**2 + D/Ei) tW = (3.0 * hb**2 * (kzi + kzf)**2 * TS) / (kb * alpha) # =2W # apply correction sig[:, 2] *= np.exp(-1*tW) # save results for peak in signals signals[index, :] = sig # g1, g1, corrected signal signals[index, 3] = peak.area # normalize by the largest peak value signals[:, 2] /= signals[:, 2].max() # simualtion signals[:, 3] /= signals[:, 3].max() # measured area return signals
[docs]def sigma_area(param): """ calculate error function between measurement and simulation """ # set nice value to 10 if not on windows try: sys.getwindowsversion() except AttributeError: isWindows = False else: isWindows = True if not isWindows: os.nice(10) # unpack parameters (spectra, h, D, xi, a, z, gclose, alpha, m, TSD, st_p_wl) = param # correct units D = D*meV xi = xi/A # list of absolute errors abs(sim-area) deltas = [] # get alpha for a given potential depth D alpha_D = get_alpha(D, *alpha) # calculate the absolute error between the simulation and the measurement # for one spectrum and write it to delta for spectrum in spectra: signal = signal_per_spectrum(spectrum, D, xi, h, a, z, gclose, alpha_D, m, TSD, st_p_wl) delta = abs(signal[:, 2]-signal[:, 3]) deltas.append(delta) # calculate least square error ? maybe not the correct term deltas = np.hstack(deltas) sigma = np.sqrt(sum(deltas**2)/(len(deltas)-len(spectra))) return sigma