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