--- /dev/null
+#!/usr/bin/python
+
+# Analyse stuff
+
+from process import *
+import odict
+import time
+
+tcsdir = "../../research/TCS/"
+imagedir = "analysis/"
+os.system("mkdir -p " + imagedir.strip("/"))
+
+defaults = [3,2,1e-8,1000,0,1]
+
+analysis = odict.odict([
+ ("sample1_focus1_si" , ["2012-08-09/074447.dat"] + defaults),
+ ("sample1_focus1_au3min", ["2012-08-09/121616.dat"] + defaults),
+ ("sample1_focus1_au8min", ["2012-08-09/161028.dat"] + defaults),
+ ("sample1_focus1_au8min_a15h", ["2012-08-10/065228.dat"]+defaults),
+ ("sample1_focus1_au18min", ["2012-08-10/104218.dat"] + defaults),
+ ("sample1_focus1_au18min_a96h", ["2012-08-14/094851.dat"] + defaults),
+ ("sample1_focus1_au33min",["2012-08-16/124728.dat"] + defaults),
+ ("sample1_focus1_au50min", ["2012-08-17/074929.dat"] + defaults),
+ ("sample1_focus1_au50min_a1h", ["2012-08-17/084544.dat"] + defaults),
+ ("sample1_focus1_au83min", ["2012-08-19/172754.dat"] + defaults),
+ ("sample1_focus1_au100min", ["2012-08-20/au_on_si(50min)/103534.dat"] + defaults),
+ ("sample1_focus1_au100min_a2h", ["2012-08-20/124236.dat"] + defaults),
+ ("sample1_focus1_au125min", ["2012-08-20/153656.dat"] + defaults),
+ ("sample1_focus1_au125min_a7d", ["2012-08-27/124551.dat"] + defaults),
+ ("sample1_focus1_au145min", ["2012-08-27/141116.dat"] + defaults),
+ ("sample1_focus2_au145min_avg", ["2012-09-07/average.dat"] + defaults),
+ ("sample1_focus2_au170min", ["2012-09-10/+30min_shinyAu/bad/"] + defaults),
+ ("sample_holder1_focus3", ["2012-09-12/sweeps/"] + defaults),
+ ("sample1_focus4_au170min", ["2012-09-13/control/"] + defaults),
+ ("sample1_focus4_blackau30min_on_au170min", ["2012-09-13/BlackAu/155950.dat"] + defaults),
+ ("sample1_focus4_blackau30min_on_au170min", ["2012-09-13/BlackAu/155434.dat"] + defaults),
+ ("sample1_focus4_blackau30min_on_au170min", ["2012-09-13/BlackAu/154918.dat"] + defaults),
+ ("sample1_focus4_au30min_on_blackau30min_on_au170min", ["2012-09-13/Au on BlackAu/"] + defaults),
+ ("sample2_focusA_si", ["2012-09-20/171922.dat"] + defaults),
+ ("sample2_focusA_si_a7d", ["2012-09-27/Si/095131.dat"] + defaults),
+ ("sample2_focusB_si_a7d", ["2012-09-27/Si (focused)/140416.dat"] + defaults),
+ ("sample2_focusB_au30min_avg", ["2012-10-01/"] + defaults),
+ ("sample2_focusB_au100min_avg", ["2012-10-03/Au (100min) on Si/"] + defaults),
+ ("sample2_focusB_au140min_avg", ["2012-10-03/Au (140min) on Si/"] + defaults),
+ ("sample2_focusB_au200min_avg", ["2012-10-03/Au (200min) on Si/"] + defaults),
+ ("sample2_focusB_au230min_avg", ["2012-10-04/Au (230min) on Si/"] + defaults),
+ ("sample2_focusC_au230min_avg", ["2012-10-04/Au (230min) on Si - Refocused/"] + defaults),
+ ("sample2_focusC_au230min_avg_a5d", ["2012-10-09/Au (230min) on Si/"] + defaults),
+ ("sample2_focusC_au260min_avg", ["2012-10-09/Au (260min) on Si/"] + defaults),
+ ("sample2_focusD_au260min_avg", ["2012-10-11/before_refocus/"] + defaults),
+ ("sample2_focusA_au260min_avg", ["2012-10-11/attempt4/"] + defaults)
+])
+
+for p in analysis.items():
+
+ print "Analysis of " + str(p[0])
+ m = p[1][1]
+ s = p[1][2]
+ min_mse = p[1][3]
+ max_fail=p[1][4]
+ max_adjust=p[1][5]
+ divide=p[1][6]
+
+ d = GetData(tcsdir + p[1][0])
+
+ gnuplot.reset()
+ gnuplot("set term postscript colour")
+ gnuplot("set output \""+str(imagedir)+str(p[0])+".eps\"")
+
+
+
+ master_title=p[0] + " - original file " + str(p[1][0])
+
+ gnuplot("set key top left")
+ ShowData(d, master_title=master_title+"\\nRaw Data", calibrate=False, normalise=False)
+
+ d = CalibrateData(d)
+ d = MaxNormalise(d)
+
+ ShowData(d, master_title=master_title+"\\nCalibrated & Normalised", calibrate=False, normalise=False)
+
+
+ gnuplot("set key top right")
+ tcs = Derivative(d)
+ ShowTCS(tcs, raw=False, master_title=master_title+"\\nTCS (of Calibrated & Normalised Data)", show_peak=False)
+ ShowTCS(tcs, raw=False, master_title=master_title+"\\nTCS with peaks (thick) and inflection points (dashed) indicated", show_peak=True)
+
+
+
+ smooth_plots = []
+ smoothed = [tcs]
+ for i in range(0, s+3):
+ smoothed.append(Smooth(smoothed[i], m=m))
+ smooth_plots += ShowTCS(smoothed[i+1], raw=False, show_peak=False, plot=None, with_="l", title="smooth x "+str(i))
+ time.sleep(0.1)
+
+ smooth_plots += ShowTCS(smoothed[0], raw=False, show_peak=False, plot=None, title="unsmoothed", with_="l lt -1 lw 1")
+
+ gnuplot("set title \""+master_title+"\\nTCS generated when "+str(m)+"pt smoothing algorithm is applied to I(E) curve\\n(calibrated and normalised)\"")
+
+ gnuplot.plot(*smooth_plots)
+
+ ShowTCS(smoothed[s], raw=False, show_peak=False, title="smooth x "+str(s), master_title=master_title+"\\nA selected level of smoothing ("+ str(s) + ")")
+
+
+ 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")
+
+ if os.path.exists(imagedir+p[0]+".fit"):
+ fit = [LoadFit(imagedir+p[0]+".fit")]
+ def temp(e):
+ total = 0.0
+ for f in fit[0]:
+ total += f[0] * gaussian(e - f[1], f[2])
+ return total
+
+ model = table(lambda e : [0.0, e, temp(e), 0.0], 0.0, 16.8, 16.8/401)
+ fit.append(model)
+ plotList = ShowTCS(smoothed[s], raw=False, plot=None)
+ plotList.append(Gnuplot.Data(model, using="2:3", with_="l lt 3"))
+ fit.append(plotList)
+ else:
+
+ #print "Not doing fit"
+ #continue
+ fit = FitTCS(smoothed[s], min_mse=min_mse, max_fail=max_fail, max_adjust=max_adjust,divide=divide, plot=None)
+ SaveFit(imagedir+p[0]+".fit", fit[0])
+
+ gnuplot("set title \""+master_title+"\\nFit gaussian peaks to smoothed TCS\"")
+ gnuplot.plot(*(fit[2]))
+
+
+ ShowTCS(fit[1], raw=False, show_peak=True, title="model", master_title=master_title+"\\nLocation of peaks in model")
+
+ components = []
+ for peak in fit[0]:
+ 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"))
+
+ gnuplot("set title \""+master_title+"\\nIndividual gaussian peaks from fit\"")
+ ymax = max(fit[0], key = lambda e : e[0] * gaussian(0, e[2]))
+ ymin = min(fit[0], key = lambda e : e[0] * gaussian(0, e[2]))
+ ymin = ymin[0] * gaussian(0, ymin[2])
+ ymax = ymax[0] * gaussian(0, ymax[2])
+
+ gnuplot("set yrange ["+str(ymin)+":"+str(ymax)+"]")
+ gnuplot.plot(*components)
+
+ elastic = fit[0][0]
+ for i in range(1, len(fit[0])):
+ peak = fit[0][i]
+ #peak[1] = elastic[1] - peak[1]
+ peak[2] = abs(peak[2]**2.0 - elastic[2]**2.0)**0.5
+
+ components = []
+ for peak in fit[0]:
+ c = table(lambda e : [0.0, e, peak[0] * gaussian(e - peak[1], peak[2]),0.0], 0, 16.8, 16.8/len(tcs))
+ components.append(Gnuplot.Data(c, using="2:3", title="", with_="l"))
+ peak = SmoothPeakFind(c, ap=DoNothing, stop=1, inflection=1)
+ components += PlotPeaks(peak,with_="l lt -1", plot=None)
+
+ gnuplot("set title \""+master_title+"\\nDeconvolved gaussian peaks from fit - Peaks shown\"")
+ ymax = max(fit[0], key = lambda e : e[0] * gaussian(0, e[2]))
+ ymin = min(fit[0], key = lambda e : e[0] * gaussian(0, e[2]))
+ ymin = ymin[0] * gaussian(0, ymin[2])
+ ymax = ymax[0] * gaussian(0, ymax[2])
+
+
+ gnuplot("set yrange ["+str(ymin)+":"+str(ymax)+"]")
+
+ gnuplot.plot(*components)
+
+ # Integrate the model
+ I = table(lambda e : [0.0, fit[1][e][1], IntegrateTCS(fit[1], 0, e),0.0], 0, len(fit[1])-1, 1)
+
+ plotItems = ShowData(I, calibrate=False, normalise=False, title="model", with_="l lt 3", plot=None)
+ plotItems += (ShowData(d, calibrate=False, normalise=False, title="data", with_="lp lt 1",plot=None))
+ time.sleep(0.1)
+ gnuplot("set yrange [0:1]")
+ gnuplot("set key top left")
+ gnuplot("set title \""+master_title+"\\nReconstructed curve from TCS model vs original data")
+ gnuplot.plot(*plotItems)
+ time.sleep(0.1)
+
+
+
git.ucc.asn.au / matches / honours.git / blobdiff
UCC git Repository :: git.ucc.asn.au