

In [65]:

    
#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
import sys
sys.path.append("../")
sys.path.append("../lmfit")
from lmfit import lmfit
import scipy.constants as const
#import the real thing
from tof_analysis import *
from helper import Moments
#some typesetting preferences
matplotlib.rcParams.update({'font.size': 14, 'font.family': 'serif'})









    



Populating the interactive namespace from numpy and matplotlib






    



WARNING: pylab import has clobbered these variables: ['helper']
`%pylab --no-import-all` prevents importing * from pylab and numpy

Determination of Surface Position

Assuming we have measured the transmitted IR Laser Power $I$ as a function of the Surface Position $Y_\mathrm{S}$.

The "SurfacePosMeas" class is designed to handle such a dataset. Its constructor's call signature is

SurfacePosMeas(Y=[], Power=[], IRPos=7)

where IRPos is the horizontal Position of the IR Laser on the Beamtool (in div) during the Surface Position Measurement.



In [66]:

    
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

As one can see, the class automatically normalizes the data and fits it to a conjugate errorfunction yielding the surface position. Called like above, the Constructor recieves an IR Laser position of +7 div on the beamtool as default value. This can be changed easily:



In [67]:

    
PosMeas2 = SurfacePosMeas(Y, Power, IRPos=6)

Of course nothing changes in the plot or the fit since IRPos does not matter here, it comes into play when calculating the different distances occuring in the actual experiment.



In [68]:

    
print PosMeas1.Y0 #The Surface Position determined from the fit
print PosMeas1.IRPos
print PosMeas1.Y









    



34.1593864228
7
[ 33.   33.2  33.4  33.6  33.8  34.   34.2  34.4  34.6  34.8  35. ]



In [69]:

    
print PosMeas1.Irel









    



[ 100.          100.           99.40828402   98.22485207   95.26627219
   76.92307692   40.82840237   14.79289941    6.62721893    3.72781065
    2.36686391]

The SurfacePosMeas class instance contains its own lmfit instance from which all attributes and properties are available:



In [70]:

    
print PosMeas1.fit.P #print the fit parameters ...
print PosMeas1.fit.StdDev #and their standard deviations









    



{'a': 3.0890480376442624, 'y0': 34.159386422793162, 'Imax': 49.963999474750516}
{'a': 7.4458529788447025, 'y0': 0.40536663333141276, 'Imax': 22.471629132866614}



In [71]:

    
PosMeas1.fit.report()









    



===========================================
Report:
-------------------------------------------
Initial set of parameters:
{'a': 1, 'y0': 35.0, 'Imax': 50}
Number of degrees of freedom (ndf): 7

Results:
Number of function evaluations: 49
      Sum of residuals (Chi^2): 56430.605614
 Variance of Chi^2 (Chi^2/ndf): 8061.515088
RMS of Chi^2 (sqrt(Chi^2/ndf)): 89.785940

Covariance Matrix:
[[  5.54407266e+01  -6.24436097e+01   5.07144220e-01]
 [ -6.24436097e+01   5.04974116e+02  -4.10072250e+00]
 [  5.07144220e-01  -4.10072250e+00   1.64322107e-01]]

Final set of parameters: 
a = 3.089048 +/- 7.445853
y0 = 34.159386 +/- 0.405367
Imax = 49.963999 +/- 22.471629
===========================================



In [72]:

    
print PosMeas1.fit.full_results
print 'RMS of Chi^2: %f' % PosMeas1.fit.full_results['RMSChi2'] # access Dictionary value for key 'RMSChi2'









    



{'VarRes': 8061.5150877529504, 'CovMatr': array([[  5.54407266e+01,  -6.24436097e+01,   5.07144220e-01],
       [ -6.24436097e+01,   5.04974116e+02,  -4.10072250e+00],
       [  5.07144220e-01,  -4.10072250e+00,   1.64322107e-01]]), 'MINPACKMsg': 'Both actual and predicted relative reductions in the sum of squares\n  are at most 0.000000', 'StdDev': {'a': 7.4458529788447025, 'y0': 0.40536663333141276, 'Imax': 22.471629132866614}, 'Parameters': {'a': 3.0890480376442624, 'y0': 34.159386422793162, 'Imax': 49.963999474750516}, 'RMSChi2': 89.78594036792704, 'Chi2': 56430.605614270651, 'fvec': array([ 0.07202147,  0.07338777, -0.47432369, -0.97683659,  1.15476748,
        1.28989313, -2.09945818,  0.14369686,  3.91678341,  3.47134631,
        2.35485286]), 'nfev': 49, 'Residuals': array([ 0.07202147,  0.07338777, -0.47432369, -0.97683659,  1.15476748,
        1.28989313, -2.09945818,  0.14369686,  3.91678341,  3.47134631,
        2.35485286])}
RMS of Chi^2: 89.785940



In [73]:

    
PosMeas1.fit.plot() #Analysis Plots for the fit (kind of silly here...)

Beamer Setup

The Positions and finally the distances in the tagging experiment are calculated by the "TaggingSetup" class. It is constructed as

TaggingSetup(self, Positions={'IR': (7,0), 'REMPI': 0,'Center ZRM': 5.1, 'ZRM': 3.5, 'Surface Y': None}, SurfacePosMeas=None, visualize=False):

The "Position" argument is a dictionary containing all relevant Positions. If "SurfacePosMeas" stays "None", then tagging of the molecules on the way back is assumed and the direct distance between IR and REMPI Laser is calculated.

If "SurfacePosMeas" argument is an instance of the SurfacePosMeas class (like the one in the upper section) then the actual Surface Position is and all other distances are calculated.



In [74]:

    
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 [75]:

    
Setup1.savefig('Setup1.pdf') # save the graph



In [76]:

    
print Setup1.FlightLength #the actual flight distance









    



15.6753326058



In [77]:

    
Setup1.show() #write everything to the standard output









    



Positions (div on beamtool):
REMPI Beam: 0.000000, 2.032000
IR Beam:    7.000000, 0.000000

Distances (mm): 
IR-S = 0.559386
IR-S-MPI = 15.675333
MPI-S = 15.115946



In [78]:

    
print Setup1.distances
print Setup1.distances['IR-S-MPI'] == Setup1.FlightLength









    



{'IR-S': 0.55938642279316042, 'IR-S-MPI': 15.675332605821659, 'MPI-S': 15.115946183028498}
True

TOF Analysis

The TOFSpectrum class

We will now look at some data from June 4th 2013

First the Setup:



In [79]:

    
Positions={'IR': (7,0), 'REMPI': -4, 'Center ZRM': 5.1, 'ZRM': 3.5} #Create dictionary with Positions
Setup = TaggingSetup(Positions, visualize=True) #Create instance of TaggingSetup class

The TOFSpectrum class is designed to handle one TOF data file and analyze it. It's constructor's signatur is

TOFSpectrum(Setup, filename, IRDelay, subtract=None, mode='Flux_vs_TOF', mass=28, verbose=False, plot=False, fit=True, baseavg=20)

verbose and plot simply trigger the the report and the plot from the lmfit class, fit argument controls whether the Data is automatically fitted or not.

Also available as modes are: 'Flux_vs_TOF', 'Flux_vs_v', 'Flux_vs_E', 'Density_vs_v'

The Data is treated as follows:

The raw data is loaded from file specified by filename
TOF is converted to $\mu$s, IRDelay is subtracted and an offset of 400 $\mu$s is added
The Baseline is subtracted from the Signal
The resulting Table is sorted along the TOF axis
If there are negative TOF values, these are omitted
TOF data gets inverted to yield the velocity $v = l/TOF$ in m/s and energy $E=\frac{1}{2}mv^2$ in eV. The length $l$ is taken from the Setup instance passed to the constructor
A baseline vertical offset is subtracted from the Signal. This offset is taken as the average of the first (baseavg=20) data points.
A Density to Flux transformation yields the (normalized) Flux $F = \mathrm{SIGNAL} \cdot v$

Fitting:

Data is fit through the Levenberg-Marquardt algorithm. The following functions are used depending on the mode:

$$F(t) = F_0 \ \left( \frac{l}{t} \right)^4\mathrm{e}^{- \left(\frac{l}{\alpha}\right)^2 \left(\frac{1}{t} - \frac{1}{t_0} \right)^2}$$$$F(v) = F_0\ v^3\ \mathrm{e}^{-(v-v_0)^2/\alpha} $$$$ F(E) = \frac{F_0}{m^2} E\ \mathrm{e}^{- \frac{2}{m\alpha^2} \left( \sqrt{E} - \sqrt{E_0} \right)}$$$$ I(v) = I_0\ v^2\ \mathrm{e}^{-(v-v_0)^2/\alpha} $$



In [80]:

    
file0406002 = TaggingTOF(Setup, '../testdata/0406002.dat', 663, plot=True, verbose=True)









    



===========================================
Report:
-------------------------------------------
Initial set of parameters:
{'F0': 1.0462185090029203, 'alpha': 7.0, 'x0': 17.0}
Number of degrees of freedom (ndf): 197

Results:
Number of function evaluations: 190
      Sum of residuals (Chi^2): 20.421849
 Variance of Chi^2 (Chi^2/ndf): 0.103664
RMS of Chi^2 (sqrt(Chi^2/ndf)): 0.321969

Covariance Matrix:
[[  1.82025176e-03  -2.37869158e-04   8.28111398e-03]
 [ -2.37869158e-04   1.32604681e-03   1.88622236e-02]
 [  8.28111398e-03   1.88622236e-02   4.55015560e-01]]

Final set of parameters: 
F0 = 0.454557 +/- 0.042664
alpha = 0.326061 +/- 0.036415
x0 = 19.891484 +/- 0.674548
===========================================

The Setup is also always available...



In [81]:

    
file0406002.Setup.visualize()

Same is true for the fit...



In [82]:

    
file0406002.Fit.P









    Out[82]:





{'F0': 0.4545570170170416,
 'alpha': 0.32606147188185797,
 'x0': 19.891483735045217}



In [83]:

    
file0406002.plot() #The plot function of TOFSpectrum class ...



In [84]:

    
file0406002.Fit.plot() #vs. the plot function of lmfit class

Same file, different mode...



In [85]:

    
file0406002 = TaggingTOF(Setup, '../testdata/0406002.dat', 663, mode='Flux_vs_E')
file0406002.plot()









    



../tof_analysis.py:483: RuntimeWarning: invalid value encountered in sqrt
  (sqrt(x)- sqrt(x0))**2 )

The TOFSpectrum class also holds an instance of the figure and axes objects from matplotlib...



In [86]:

    
file0406002.savefig('file0406002.svg') #save graph as svg vector image
file0406002.savefig('file0406002.pdf') # or pdf ... or eps... or... c.f. matplotlib reference for available graphic formats



In [87]:

    
file0406002.set_mode('Flux_vs_v') #change mode (but its better to create a new TOFSpectrum instance for that)
file0406002.fit()
file0406002.plot()



In [88]:

    
file0406002.set_mode('Flux_vs_E') #switch back to Energy
file0406002.fit()
file0406002.plot() #instances of matplotlib.figure and axes exist only when TOFSpectrum.plot() was calles before at some point
file0406002.axes.set_xlim((0,3)) #set xaxis limits. TOFSpectrum.axes gives you full access to the graphs properties
display(file0406002.fig) # show figure object in inline mode

The first three moments around 0 of the fit function can be calculated. Those are derived by numerical integration from 0 to $\infty$ using scipy.integrate.quad(). $$ \langle x^n \rangle = \int\limits_0^\infty x^n f(x) \ \mathrm{d} x $$ where $n=0,1,2$.



In [89]:

    
moments, mean, variance = file0406002.Moments
print moments # moments with numerical integration errors









    



0.1395912012

The Mean and the Mean Squared values of the distributions can be calculated:

$$\langle x \rangle = \frac{\int\limits_0^\infty xf(x) \ \mathrm{d} x}{\int\limits_0^\infty f(x) \ \mathrm{d}x} $$$$\langle x^2 \rangle = \frac{\int\limits_0^\infty x^2f(x) \ \mathrm{d} x}{\int\limits_0^\infty f(x) \ \mathrm{d}x} $$



In [90]:

    
print mean    
print variance
RMS = sqrt(variance)
print RMS









    



0.0255456687604
0.00503911898109
0.070986752152

This works as well in other modes...



In [91]:

    
file0406002_v = TaggingTOF(Setup, '../testdata/0406002.dat', 663, mode='Flux_vs_v')
file0406002_v.plot()



In [92]:

    
print file0406002_v.Moments[1]
print sqrt(file0406002_v.Moments[2])









    



686616.022305
29407.743744



In [93]:

    
file0406002_E = TaggingTOF(Setup, '../testdata/0406002.dat', 663, mode='Flux_vs_E')
file0406002_E.plot()
file0406002_E.axes.set_xlim((0,1.5)) #set xaxis limits. TOFSpectrum.axes gives you full access to the graphs properties
display(file0406002_E.fig) # show figure object in inline mode



In [94]:

    
print file0406002_E.Fit.P
print file0406002_E.Fit.StdDev









    



{'F0': 4118.3739127393401, 'alpha': -0.035997925870835799, 'x0': 0.16179907382728695}
{'F0': 359.77440999530597, 'alpha': 0.0038390265168828444, 'x0': 0.0088232391577434018}



In [95]:

    
P_Mean, P_StdDev, outlist = file0406002_E.Fit.bootstrap(n=500, plot=True)
print P_Mean
print P_StdDev
fitfunc = file0406002.get_fitfunc() #get the current fitfunction
E_Mean = np.array([None for i in outlist])
E_RMS = np.array([None for i in outlist])
i = 0
for finalParamsDict in outlist: # iterate over the final parameters of the bootstrap fits
    zeroth, first, second = Moments(fitfunc, [0,1,2], args=finalParamsDict)
    E_Mean[i] = first[0]/zeroth[0]
    E_RMS[i] = second[0]/zeroth[0]
    i += 1
print 'Mean Energy: (%f +/- %f) eV' %(E_Mean.mean(), E_Mean.std()) # maybe you want to use the standard error of the mean here instead of the standard deviation
print 'RMS Energy: (%f +/- %f) eV' %(E_RMS.mean(), E_RMS.std())









    



Bootstrapping:






    












    



{'F0': 4040.3626668556744, 'alpha': -0.034586761737545754, 'x0': 0.16436420026225412}
{'F0': 69.681639165076291, 'alpha': 0.0006761963846957271, 'x0': 0.0015430194479016627}
Mean Energy: (0.220778 +/- 0.002296) eV
RMS Energy: (0.055234 +/- 0.001264) eV



In [96]:

    
file0406003 = TaggingTOF(Setup, '../testdata/0406003.dat', 663)
file0406003.plot()

Using TOFSpectrum for manual analysis:



In [97]:

    
TOF, Flux = file0406003.get_values()
TOFslice, Fluxslice = file0406003.get_values(limits=(10,30))
plot(TOF, Flux, 'bo')
plot(TOFslice, Fluxslice, 'ro')
fitfunc = file0406003.get_fitfunc()
print fitfunc.func_code.co_varnames #inspect live object to see the variable names









    



('self', 'x', 'F0', 'alpha', 'x0')



In [98]:

    
file0406003fit = lmfit(fitfunc, TOFslice, Fluxslice, p0={'F0':1, 'alpha': -0.2, 'x0':20}, verbose=False, plot=True, plot_options={'acf':False, 'lagplot':False})



In [99]:

    
zeroth, first, second = Moments(file0406003fit, [0,1,2])
mean = first[0]/zeroth[0]



In [100]:

    
fig1 = figure(1)
plot(TOF, Flux, 'bo')
ylim((-0.01, 1.05))
xlim((0,50))
xlabel(r'$\mathrm{TOF}\ / \ \mu\mathrm{s}$')
ylabel(r'$\mathrm{Normalized\ Flux}$')
plot(TOFslice, Fluxslice, 'ro')
x = linspace(TOF[0], TOF[-1], 1000)
plot(x, file0406003fit(x), 'r-', label=r'Fit')
legend(loc='best')
bar(mean, file0406003fit(mean), width=0.01)
string = r'$\langle \mathrm{TOF} \rangle =\ %2.2f\ \mu\mathrm{s}$' % mean
annotate(string, xy=(25,0.6), xycoords='data')
savefig('0406003manual.pdf')
plt.show()

Dealing with Background:



In [101]:

    
Positions={'IR': (7,0), 'REMPI': -4, 'Center ZRM': 5.1, 'ZRM': 3.5, 'Surface Y': 33.6} #Create dictionary with Positions
Setup3 = TaggingSetup(Positions, visualize=True) #Create instance of TaggingSetup class
file0606028 = TaggingTOF(Setup, '../testdata/0606028.dat', 654, fit=False)
file0606029 = TaggingTOF(Setup, '../testdata/0606029.dat', 654, fit=False)



In [102]:

    
file0606028.RawData.plot()



In [103]:

    
file0606029.RawData.plot()



In [104]:

    
TOF = file0606029.RawData.TOF
SIG = file0606029.RawData.Baseline-file0606029.RawData.Signal
bgfit = polyfit(TOF, SIG, 5)
print bgfit
plot(TOF, SIG)
x = linspace(xlim()[0], xlim()[1], 1000)
plot(x, polyval(bgfit, x))
show()









    



[ -1.49151672e-27   2.04953191e-21  -9.53294332e-16   1.51653693e-10
   2.45522183e-06  -1.31037752e+00]



In [105]:

    
Sig = file0606028.RawData.Baseline-file0606028.RawData.Signal - polyval(bgfit, file0606029.RawData.TOF)
plot(file0606028.RawData.TOF, Sig, 'k-')
show()



In [106]:

    
savetxt('temp.dat', transpose(array([file0606028.RawData.TOF, -Sig, zeros(len(Sig))]))) # save to file, third column with zeros must be appended



In [107]:

    
file0606028corr = TaggingTOF(Setup3, 'temp.dat', 654, fit=False)



In [108]:

    
TOF, Flux = file0606028corr.TOF
TOFslice, Fluxslice = file0606028corr.get_values(limits=(10,160))
plot(TOF, Flux, 'bo')
plot(TOFslice, Fluxslice, 'ro')
show()



In [109]:

    
fitfunc = file0606028.get_fitfunc()
file0606028corr_fit = lmfit(fitfunc, TOFslice, Fluxslice, p0={'F0':1, 'alpha':-0.2, 'x0':20}, verbose=False, plot=True)



In [110]:

    
file0606028corr2 = TaggingTOF(Setup3, 'temp.dat', 654, fit=False)
file0606028corr2.fit(p0={'F0':1, 'alpha':-0.2, 'x0':20}, plot=False)
file0606028corr2.plot()

The MultiTOF class



In [111]:

    
files = [u'../testdata/0406012.dat',\
      u'../testdata/0406013.dat',\
      u'../testdata/0406014.dat',\
      u'../testdata/0406015.dat',\
      u'../testdata/0406016.dat',\
      u'../testdata/0406017.dat',\
      u'../testdata/0406018.dat',\
      u'../testdata/0406019.dat',\
      u'../testdata/0406020.dat',\
      u'../testdata/0406021.dat',\
      u'../testdata/0406022.dat',\
      u'../testdata/0406023.dat']
descriptions = [r'$100\ ^\circ \mathrm{C}$',\
      r'$200\ ^\circ \mathrm{C}$',\
      r'$400\ ^\circ \mathrm{C}$',\
      r'$600\ ^\circ \mathrm{C}$',\
      r'$700\ ^\circ \mathrm{C}$',\
      r'$500\ ^\circ \mathrm{C}$',\
      r'$300\ ^\circ \mathrm{C}$',\
      r'$-180\ ^\circ \mathrm{C}$',\
      r'$-100\ ^\circ \mathrm{C}$',\
      r'$-50\ ^\circ \mathrm{C}$',\
      r'$0\ ^\circ \mathrm{C}$',\
      r'$30\ ^\circ \mathrm{C}$']
R6Tagging = MultiTOF(files, descriptions, Setup, 663, mode='Flux_vs_E')



In [112]:

    
R6Tagging.plot(xlim=(0,3))



In [113]:

    
R6_300deg = R6Tagging[3]
R6_300deg.plot()



In [114]:

    
R6_300deg.Fit.plot()



In [115]:

    
moments, mean, variance = R6_300deg.Moments
print mean, sqrt(variance)









    



0.0126313151914 0.0426344279872

Usefull code snippets

File opening dialog using Tkinter



In [116]:

    
from helper import SelectFileDialog

Uncomment the next Cell and run it...



In [117]:

    
#files = SelectFileDialog()
#print files

Construct filenames



In [118]:

    
from helper import construct_filename



In [119]:

    
print construct_filename('1607', 1)









    



1607001.dat



In [120]:

    
print construct_filename('1607', range(10))









    



['1607000.dat', '1607001.dat', '1607002.dat', '1607003.dat', '1607004.dat', '1607005.dat', '1607006.dat', '1607007.dat', '1607008.dat', '1607009.dat']

Clean this up



In [121]:

    
ls









    



0406003manual.pdf    file0406002.svg  simulations.ipynb
angular_distr.ipynb  newtof.ipynb     temp.dat
file0406002.pdf      Setup1.pdf       tof_analysis.ipynb



In [122]:

    
rm -f file0406002.pdf



In [123]:

    
rm -f file0406002.svg



In [124]:

    
rm -f Setup1.pdf



In [125]:

    
rm -f temp.dat



In [126]:

    
rm -f 0406003manual.pdf



In [126]: