9 tcsdir = "../../research/TCS/"
10 imagedir = "analysis/"
11 os.system("mkdir -p " + imagedir.strip("/"))
13 defaults = [3,2,1e-8,1000,0,1]
15 analysis = odict.odict([
16 ("sample1_focus1_si" , ["2012-08-09/074447.dat"] + defaults),
17 ("sample1_focus1_au3min", ["2012-08-09/121616.dat"] + defaults),
18 ("sample1_focus1_au8min", ["2012-08-09/161028.dat"] + defaults),
19 ("sample1_focus1_au8min_a15h", ["2012-08-10/065228.dat"]+defaults),
20 ("sample1_focus1_au18min", ["2012-08-10/104218.dat"] + defaults),
21 ("sample1_focus1_au18min_a96h", ["2012-08-14/094851.dat"] + defaults),
22 ("sample1_focus1_au33min",["2012-08-16/124728.dat"] + defaults),
23 ("sample1_focus1_au50min", ["2012-08-17/074929.dat"] + defaults),
24 ("sample1_focus1_au50min_a1h", ["2012-08-17/084544.dat"] + defaults),
25 ("sample1_focus1_au83min", ["2012-08-19/172754.dat"] + defaults),
26 ("sample1_focus1_au100min", ["2012-08-20/au_on_si(50min)/103534.dat"] + defaults),
27 ("sample1_focus1_au100min_a2h", ["2012-08-20/124236.dat"] + defaults),
28 ("sample1_focus1_au125min", ["2012-08-20/153656.dat"] + defaults),
29 ("sample1_focus1_au125min_a7d", ["2012-08-27/124551.dat"] + defaults),
30 ("sample1_focus1_au145min", ["2012-08-27/141116.dat"] + defaults),
31 ("sample1_focus2_au145min_avg", ["2012-09-07/average.dat"] + defaults),
32 ("sample1_focus2_au170min", ["2012-09-10/+30min_shinyAu/bad/"] + defaults),
33 ("sample_holder1_focus3", ["2012-09-12/sweeps/"] + defaults),
34 ("sample1_focus4_au170min", ["2012-09-13/control/"] + defaults),
35 ("sample1_focus4_blackau30min_on_au170min", ["2012-09-13/BlackAu/155950.dat"] + defaults),
36 ("sample1_focus4_blackau30min_on_au170min", ["2012-09-13/BlackAu/155434.dat"] + defaults),
37 ("sample1_focus4_blackau30min_on_au170min", ["2012-09-13/BlackAu/154918.dat"] + defaults),
38 ("sample1_focus4_au30min_on_blackau30min_on_au170min", ["2012-09-13/Au on BlackAu/"] + defaults),
39 ("sample2_focusA_si", ["2012-09-20/171922.dat"] + defaults),
40 ("sample2_focusA_si_a7d", ["2012-09-27/Si/095131.dat"] + defaults),
41 ("sample2_focusB_si_a7d", ["2012-09-27/Si (focused)/140416.dat"] + defaults),
42 ("sample2_focusB_au30min_avg", ["2012-10-01/"] + defaults),
43 ("sample2_focusB_au100min_avg", ["2012-10-03/Au (100min) on Si/"] + defaults),
44 ("sample2_focusB_au140min_avg", ["2012-10-03/Au (140min) on Si/"] + defaults),
45 ("sample2_focusB_au200min_avg", ["2012-10-03/Au (200min) on Si/"] + defaults),
46 ("sample2_focusB_au230min_avg", ["2012-10-04/Au (230min) on Si/"] + defaults),
47 ("sample2_focusC_au230min_avg", ["2012-10-04/Au (230min) on Si - Refocused/"] + defaults),
48 ("sample2_focusC_au230min_avg_a5d", ["2012-10-09/Au (230min) on Si/"] + defaults),
49 ("sample2_focusC_au260min_avg", ["2012-10-09/Au (260min) on Si/"] + defaults),
50 ("sample2_focusD_au260min_avg", ["2012-10-11/before_refocus/"] + defaults),
51 ("sample2_focusA_au260min_avg", ["2012-10-11/attempt4/"] + defaults)
54 for p in analysis.items():
56 print "Analysis of " + str(p[0])
64 d = GetData(tcsdir + p[1][0])
67 gnuplot("set term postscript colour")
68 gnuplot("set output \""+str(imagedir)+str(p[0])+".eps\"")
72 master_title=p[0] + " - original file " + str(p[1][0])
74 gnuplot("set key top left")
75 ShowData(d, master_title=master_title+"\\nRaw Data", calibrate=False, normalise=False)
80 ShowData(d, master_title=master_title+"\\nCalibrated & Normalised", calibrate=False, normalise=False)
83 gnuplot("set key top right")
85 ShowTCS(tcs, raw=False, master_title=master_title+"\\nTCS (of Calibrated & Normalised Data)", show_peak=False)
86 ShowTCS(tcs, raw=False, master_title=master_title+"\\nTCS with peaks (thick) and inflection points (dashed) indicated", show_peak=True)
92 for i in range(0, s+3):
93 smoothed.append(Smooth(smoothed[i], m=m))
94 smooth_plots += ShowTCS(smoothed[i+1], raw=False, show_peak=False, plot=None, with_="l", title="smooth x "+str(i))
97 smooth_plots += ShowTCS(smoothed[0], raw=False, show_peak=False, plot=None, title="unsmoothed", with_="l lt -1 lw 1")
99 gnuplot("set title \""+master_title+"\\nTCS generated when "+str(m)+"pt smoothing algorithm is applied to I(E) curve\\n(calibrated and normalised)\"")
101 gnuplot.plot(*smooth_plots)
103 ShowTCS(smoothed[s], raw=False, show_peak=False, title="smooth x "+str(s), master_title=master_title+"\\nA selected level of smoothing ("+ str(s) + ")")
106 ShowTCS(smoothed[s], raw=False, show_peak=True, title="smooth x "+str(s), master_title=master_title+"\\nPeaks and inflection points of TCS from selected smoothing level")
108 if os.path.exists(imagedir+p[0]+".fit"):
109 fit = [LoadFit(imagedir+p[0]+".fit")]
113 total += f[0] * gaussian(e - f[1], f[2])
116 model = table(lambda e : [0.0, e, temp(e), 0.0], 0.0, 16.8, 16.8/401)
118 plotList = ShowTCS(smoothed[s], raw=False, plot=None)
119 plotList.append(Gnuplot.Data(model, using="2:3", with_="l lt 3"))
123 #print "Not doing fit"
125 fit = FitTCS(smoothed[s], min_mse=min_mse, max_fail=max_fail, max_adjust=max_adjust,divide=divide, plot=None)
126 SaveFit(imagedir+p[0]+".fit", fit[0])
128 gnuplot("set title \""+master_title+"\\nFit gaussian peaks to smoothed TCS\"")
129 gnuplot.plot(*(fit[2]))
132 ShowTCS(fit[1], raw=False, show_peak=True, title="model", master_title=master_title+"\\nLocation of peaks in model")
136 components.append(Gnuplot.Data(table(lambda e : [0.0, e, peak[0] * gaussian(e - peak[1], peak[2]), 0.0], 0, 16.8, 16.8/len(tcs)), using="2:3", title="", with_="l"))
138 gnuplot("set title \""+master_title+"\\nIndividual gaussian peaks from fit\"")
139 ymax = max(fit[0], key = lambda e : e[0] * gaussian(0, e[2]))
140 ymin = min(fit[0], key = lambda e : e[0] * gaussian(0, e[2]))
141 ymin = ymin[0] * gaussian(0, ymin[2])
142 ymax = ymax[0] * gaussian(0, ymax[2])
144 gnuplot("set yrange ["+str(ymin)+":"+str(ymax)+"]")
145 gnuplot.plot(*components)
148 for i in range(1, len(fit[0])):
150 #peak[1] = elastic[1] - peak[1]
151 peak[2] = abs(peak[2]**2.0 - elastic[2]**2.0)**0.5
155 c = table(lambda e : [0.0, e, peak[0] * gaussian(e - peak[1], peak[2]),0.0], 0, 16.8, 16.8/len(tcs))
156 components.append(Gnuplot.Data(c, using="2:3", title="", with_="l"))
157 peak = SmoothPeakFind(c, ap=DoNothing, stop=1, inflection=1)
158 components += PlotPeaks(peak,with_="l lt -1", plot=None)
160 gnuplot("set title \""+master_title+"\\nDeconvolved gaussian peaks from fit - Peaks shown\"")
161 ymax = max(fit[0], key = lambda e : e[0] * gaussian(0, e[2]))
162 ymin = min(fit[0], key = lambda e : e[0] * gaussian(0, e[2]))
163 ymin = ymin[0] * gaussian(0, ymin[2])
164 ymax = ymax[0] * gaussian(0, ymax[2])
167 gnuplot("set yrange ["+str(ymin)+":"+str(ymax)+"]")
169 gnuplot.plot(*components)
171 # Integrate the model
172 I = table(lambda e : [0.0, fit[1][e][1], IntegrateTCS(fit[1], 0, e),0.0], 0, len(fit[1])-1, 1)
174 plotItems = ShowData(I, calibrate=False, normalise=False, title="model", with_="l lt 3", plot=None)
175 plotItems += (ShowData(d, calibrate=False, normalise=False, title="data", with_="lp lt 1",plot=None))
177 gnuplot("set yrange [0:1]")
178 gnuplot("set key top left")
179 gnuplot("set title \""+master_title+"\\nReconstructed curve from TCS model vs original data")
180 gnuplot.plot(*plotItems)