#!/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/155434.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)