In [1]:
#pylab mode
#remove the option inline to open plot windows using the selected gui backend. Induvidual graphs can stille be embedded using display(youfigure)
%pylab inline
#some imports for convenience
from __future__ import division
from lmfit import lmfit
import scipy.constants as const
#import the real thing
from tof_analysis import *
from angular_distr import REMPISetup
#some typesetting preferences
matplotlib.rcParams.update({'font.size': 14, 'font.family': 'serif'})
In [2]:
from __future__ import division
from matplotlib.pyplot import plot
import matplotlib.pyplot as plt
import numpy as np
from numpy import sqrt, cos, sin, log, exp, linspace
from matplotlib.ticker import MultipleLocator, NullFormatter
from lmfit import lmfit
from matplotlib.path import Path
import matplotlib.patches as patches
import matplotlib.gridspec as gs
import scipy.constants as const
from scipy.special import erfc
from scipy.integrate import quad
from scipy.optimize import leastsq
import helper
In [3]:
class TaggingTOF(object):
"""
...
"""
# CONSTRUCTOR
def __init__(self, Setup, filename, IRDelay,
mode='Flux_vs_TOF', mass=28,
verbose=False, plot=True, fit=True, numbasecorr=20,
Normalize=True):
"""
...
"""
self.__mass = mass
self.__filename = filename
self.__IRDelay = IRDelay
self.__setup = Setup
self.__l = Setup.FlightLength
self.__offset = 400
self.__Normalize = Normalize
self.__numbasecorr = numbasecorr
self.__RawData = RawTOFData(filename)
TOF, Signal = self.__treat_data(self.__RawData)
v, E = self.__invert(TOF)
Flux = self.__density_to_flux(Signal, v)
if Normalize:
Flux = [helper.normalize(flux, 5) for flux in Flux]
self.__TOF = (TOF, Flux[0])
self.__v = (v, Flux[1])
self.__E = (E, Flux[2])
self.set_mode(mode)
self.__moments = None
if fit:
self.fit(verbose=verbose)
if plot:
self.plot()
# PRIVATE ATTRIBUTES
def __treat_data(self, RawData):
idx = RawData.TOF.argsort()
TOF = RawData.TOF[idx]/1000 -self.__IRDelay + self.__offset
Signal = (RawData.Baseline - RawData.Signal)[idx]
Signal = helper.subtract_baseline(Signal, left=self.__numbasecorr)
idx, = np.where(TOF > 0)
numomitt = len(TOF) - len(idx)
if numomitt > 0:
print 'Warning: Omitting the first %d datapoints because TOF <= 0!' % numomitt
return TOF[idx], Signal[idx]
def __invert(self, TOF):
v = self.__l*1000 / TOF
E = 0.5 * self.__mass * v**2 * const.u/const.eV
return v, E
def __density_to_flux(self, Signal, v):
TOFFlux = Signal*v
vFlux = Signal*self.__l/v
EFlux = vFlux/(v*self.__mass)
return TOFFlux, vFlux, EFlux
def __FluxTOFfit(self, x, F0, alpha, x0):
return F0 * (self.__l/x)**4 * exp(- (self.__l/alpha)**2 *
(1/x - 1/x0)**2)
def __FluxVfit(self, x, F0, alpha, x0):
return F0 * x**3 * exp(-( (x-x0) /alpha)**2)
def __FluxEfit(self, x, F0, alpha, x0):
return F0 * x/self.__mass**2 * exp(- 2/(self.__mass*alpha**2) *
(sqrt(x)- sqrt(x0))**2 )
def __DensityVfit(self, x, F0, alpha, x0):
return F0 * x**2 * exp(- ((x-x0)/alpha)**2)
def __estimate_p0(self, x, y):
maxpos = np.where(y == y.max())[0][0]
x0 = x[maxpos]
F0 = y[maxpos]
i = maxpos
while y[i] > F0/2:
i += 1
right = i
i = maxpos
while y[i] > F0/2:
i -= 1
left = i
alpha = x[right] - x[left]
return {'F0':F0, 'x0':x0, 'alpha':alpha}
def __get_moments(self):
if self.__moments == None:
zeroth, first, second = helper.Moments(self.__fit, n=[0,1,2], limits=(0, np.inf))
self.__moments = ([zeroth[0], first[0], second[0]], first[0]/zeroth[0], second[0]/zeroth[0])
return self.__moments
def __call__(self, x):
return self.__fit(x)
# PUBLIC ATTRIBUTES
@property
def Mass(self):
return self.__mass
@property
def RawData(self):
return self.__RawData
@property
def TOF(self):
return self.__TOF
@property
def v(self):
return self.__v
@property
def E(self):
return self.__E
@property
def Fit(self):
return self.__fit
@property
def Setup(self):
return self.__setup
@property
def Mean(self):
return self.__get_moments()[1]
@property
def StdDev(self):
return self.__get_moments()[2]
@property
def Moments(self):
return self.__get_moments()[0]
def get_fitfunc(self):
return self.__fitargs['func']
def set_mode(self, mode):
pltset = {
'xlabel': r'$\mathrm{Time\ of\ Flight\, /\, \mu s }$',
'ylabel': r'',
}
if self.__Normalize:
pltset['ylabel'] = r'$\mathrm{Normalized}$'
if mode == 'Flux_vs_TOF':
self.__fitargs = {'func': self.__FluxTOFfit, 'data': self.__TOF}
pltset['ylabel'] = pltset['ylabel'] + r' $\mathrm{Flux}$'
elif mode == 'Flux_vs_v':
self.__fitargs = {'func': self.__FluxVfit, 'data': self.__v}
pltset['ylabel'] = pltset['ylabel'] + r' $\mathrm{Flux}$'
pltset['xlabel'] = r'$v \, / \, \mathrm{m}\, \mathrm{s}^{-1}$'
elif mode == 'Flux_vs_E':
self.__fitargs = {'func': self.__FluxEfit, 'data': self.__E}
pltset['ylabel'] = pltset['ylabel'] + r' $\mathrm{Flux}$'
pltset['xlabel'] = r'$E \, / \, \mathrm{eV}$'
else:
raise Exception('Error: Mode %s unknown!' % mode)
return None
self._plotsettings = pltset
def fit(self, func=None, p0={}, verbose=False, plot=False):
self.__moments = None
if func == None:
func = self.__fitargs['func']
x, y = self.__fitargs['data']
if p0 == {}:
p0 = self.__estimate_p0(x, y)
self.__fit = lmfit(func, x, y, p0, plot=plot, verbose=verbose)
def plot(self):
fig, axes = plt.subplots(figsize=(8,5), dpi=100)
axes.set_title(r'$\mathrm{%s}$' % self.__filename)
x, y = self.__fitargs['data']
pltset = self._plotsettings
xlim = ( 0.928*x[0], 1.02*x[-1] )
if xlim[0] > xlim[1]:
xlim = (xlim[1], xlim[0])
axes.plot(x, y, 'bo')
xspace = linspace(xlim[0], xlim[1], 1000)
axes.set_xlim(xlim)
if self.__fit:
axes.plot(xspace, self.Fit(xspace), 'r-')
axes.set_ylabel(pltset['ylabel'])
axes.set_xlabel(pltset['xlabel'])
self.fig = fig
self.axes = axes
self.savefig = self.fig.savefig
plt.show(fig)
In [4]:
Y = np.linspace(33, 35, 11) #Create a list of 11 equally spaced numbers from 33 to 35.
Power = [16.9, 16.9, 16.8, 16.6, 16.1, 13, 6.9, 2.5, 1.12, 0.63, 0.4] #Create a list of the measured IR Power corresponding to Y
PosMeas1 = SurfacePosMeas(Y, Power) #Create an instance of the SurfacePosMeas class
In [5]:
Positions={'IR': (7,0), 'REMPI': 0, 'Center ZRM': 5.1, 'ZRM': 3.5, 'Surface Y': 33.6} #Create dictionary with Positions
Setup1 = TaggingSetup(Positions, SurfacePosMeas=PosMeas1, visualize=True) #Create instance of TaggingSetup class
In [6]:
testtof = TaggingTOF(Setup1, 'testdata/1308033.dat', 606, verbose=True, mode='Flux_vs_E')
In [7]:
print testtof.Mean, testtof.StdDev, testtof.Moments
In [7]:
In [7]: