research/analysis and stuff/tcs scripts/analysis.py

   1 #!/usr/bin/python
   2
   3 # Analyse stuff
   4
   5 from process import *
   6 import odict
   7 import time
   8
   9 tcsdir = "../../research/TCS/"
  10 imagedir = "analysis/"
  11 os.system("mkdir -p " + imagedir.strip("/"))
  12
  13 defaults = [3,2,1e-8,1000,0,1]
  14
  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)
  52 ])
  53
  54 for p in analysis.items():
  55
  56         print "Analysis of " + str(p[0])
  57         m = p[1][1]
  58         s = p[1][2]
  59         min_mse = p[1][3]
  60         max_fail=p[1][4]
  61         max_adjust=p[1][5]
  62         divide=p[1][6]
  63
  64         d = GetData(tcsdir + p[1][0])
  65
  66         gnuplot.reset()
  67         gnuplot("set term postscript colour")
  68         gnuplot("set output \""+str(imagedir)+str(p[0])+".eps\"")
  69
  70
  71
  72         master_title=p[0] + " - original file " + str(p[1][0])
  73
  74         gnuplot("set key top left")
  75         ShowData(d, master_title=master_title+"\\nRaw Data", calibrate=False, normalise=False)
  76
  77         d = CalibrateData(d)
  78         d = MaxNormalise(d)
  79
  80         ShowData(d, master_title=master_title+"\\nCalibrated & Normalised", calibrate=False, normalise=False)
  81
  82
  83         gnuplot("set key top right")
  84         tcs = Derivative(d)
  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)
  87
  88
  89
  90         smooth_plots = []
  91         smoothed = [tcs]
  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))
  95                 time.sleep(0.1)
  96
  97         smooth_plots += ShowTCS(smoothed[0], raw=False, show_peak=False, plot=None, title="unsmoothed", with_="l lt -1 lw 1")
  98
  99         gnuplot("set title \""+master_title+"\\nTCS generated when "+str(m)+"pt smoothing algorithm is applied to I(E) curve\\n(calibrated and normalised)\"")
 100
 101         gnuplot.plot(*smooth_plots)
 102
 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) + ")")
 104
 105
 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")
 107
 108         if os.path.exists(imagedir+p[0]+".fit"):
 109                 fit = [LoadFit(imagedir+p[0]+".fit")]
 110                 def temp(e):
 111                         total = 0.0
 112                         for f in fit[0]:
 113                                 total += f[0] * gaussian(e - f[1], f[2])
 114                         return total
 115
 116                 model = table(lambda e : [0.0, e, temp(e), 0.0], 0.0, 16.8, 16.8/401)
 117                 fit.append(model)
 118                 plotList = ShowTCS(smoothed[s], raw=False, plot=None)
 119                 plotList.append(Gnuplot.Data(model, using="2:3", with_="l lt 3"))
 120                 fit.append(plotList)
 121         else:
 122
 123                 #print "Not doing fit"
 124                 #continue
 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])
 127
 128         gnuplot("set title \""+master_title+"\\nFit gaussian peaks to smoothed TCS\"")
 129         gnuplot.plot(*(fit[2]))
 130
 131
 132         ShowTCS(fit[1], raw=False, show_peak=True, title="model", master_title=master_title+"\\nLocation of peaks in model")
 133
 134         components = []
 135         for peak in fit[0]:
 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"))
 137
 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])
 143
 144         gnuplot("set yrange ["+str(ymin)+":"+str(ymax)+"]")
 145         gnuplot.plot(*components)
 146
 147         elastic = fit[0][0]
 148         for i in range(1, len(fit[0])):
 149                 peak = fit[0][i]
 150                 #peak[1] = elastic[1] - peak[1]
 151                 peak[2] = abs(peak[2]**2.0 - elastic[2]**2.0)**0.5
 152
 153         components = []
 154         for peak in fit[0]:
 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)
 159
 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])
 165
 166
 167         gnuplot("set yrange ["+str(ymin)+":"+str(ymax)+"]")
 168
 169         gnuplot.plot(*components)
 170
 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)
 173
 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))
 176         time.sleep(0.1)
 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)
 181         time.sleep(0.1)
 182
 183
 184