5 # @purpose Process TCS data
13 import re # Regular expressions - for removing comments
14 import odict #ordered dictionary
17 import Gnuplot, Gnuplot.funcutils
25 gnuplot = Gnuplot.Gnuplot()
28 gnuplot = Gnuplot.Gnuplot()
30 def FindDataFiles(directory=".", depth=1, result=None):
34 for f in os.listdir(directory):
35 if os.path.isdir(directory+"/"+str(f)):
37 result += FindDataFiles(directory+"/"+str(f), depth-1, result)
40 if (len(s) == 2 and s[1] == "dat"):
41 result.append(directory+"/"+str(f))
50 def DirectoryName(f, start=0,back=1):
52 return string.join(a[start:(len(a)-back)], "/")
54 def GetData(filename, key=1):
56 if type(filename) != type(""):
57 if type(filename) == type([]):
63 if os.path.isdir(filename):
64 if os.path.exists(filename.strip("/")+"/average.dat"):
65 os.remove(filename.strip("/")+"/average.dat")
66 AverageAllData(filename)
67 return GetData(filename.strip("/")+"/average.dat")
69 input_file = open(filename, "r")
71 for line in input_file:
72 line = re.sub("#.*", "", line).strip("\r\n ")
76 line = map(lambda e : float(e), line.split("\t"))
79 for i in range(len(line)):
80 data[line[key]][0][i] += line[i]
81 data[line[key]][1] += 1
83 data.update({line[key] : [line, 1]})
85 d = map(lambda e : map(lambda f : float(f) / float(e[1][1]), e[1][0]), data.items())
86 d.sort(key = lambda e : e[key])
95 def GetDataSets(directory="."):
97 for f in os.listdir(directory):
98 if os.path.isdir(directory+"/"+str(f)) == False:
99 if (len(f.split(".")) > 1 and f.split(".")[1] == "dat"):
100 d = GetData(directory+"/"+str(f))
107 def Derivative(data, a=1, b=2, sigma=None,step=1):
113 for i in range(0, len(data),step):
114 result[len(result)-1] = [d for d in data[i]]
116 dE[0] = data[i][a] - data[i-step][a]
117 dI[0] = data[i][b] - data[i-step][b]
121 if (i < len(data)-step):
122 dE[1] = data[i+step][a] - data[i][a]
123 dI[1] = data[i+step][b] - data[i][b]
128 #print str(data[i]) + " ["+str(sigma)+"] = " + str(data[i][int(abs(sigma))])
130 if dI[0] != None: dI[0] -= 0.5*data[i][int(abs(sigma))]
131 if dI[1] != None: dI[1] -= 0.5*data[i][int(abs(sigma))]
133 if dI[0] != None: dI[0] += 0.5*data[i][int(abs(sigma))]
134 if dI[1] != None: dI[1] += 0.5*data[i][int(abs(sigma))]
149 deltaI /= float(count)
150 deltaE /= float(count)
154 result[len(result)-1][b] = (deltaI / deltaE)
156 result[len(result)-1][b] = 0.0
158 result[len(result)-1][b] = 0.0
161 return result[0:len(result)-1]
163 def MaxNormalise(data, u=2):
164 result = copy.deepcopy(data)
167 maxval = max(data, key = lambda e : e[u])[u]
177 def Average(data_sets, u=1):
182 #print "Already have " + str(p[u])
184 for i in range(0, len(p)):
185 avg[p[u]][0][i] += p[i]
187 #print "Create key for " + str(p[u])
188 avg.update({p[u] : [p, 1]})
191 for i in range(0, len(avg[a][0])):
192 avg[a][0][i] /= float(avg[a][1])
193 return map(lambda e : e[1][0], sorted(avg.items(), key = lambda e : e[0]))
195 def FullWidthAtHalfMax(data, u=1):
196 maxval = max(data, key = lambda e : e[u])
197 peak = data.index(maxval)
201 for i in range(1, len(data)/2):
203 if (peak-i > 0 and data[peak-i] < 0.50*maxval):
204 lhs = data[peak-i][u]
206 if (peak+i < len(data) and data[peak+i] < 0.50*maxval):
208 if lhs != None and rhs != None:
210 if rhs == None or lhs == None:
211 return abs(data[len(data)-1][0] - data[0][0])
213 return abs(rhs - lhs)
215 def SaveData(filename, data):
216 out = open(filename, "w", 0)
218 for i in range(0, len(a)):
224 def AverageAllData(directory, save=None, normalise=True):
226 if save == None: save = directory+"/average.dat"
227 for f in FindDataFiles(directory):
233 a = Average(data_sets)
237 def CalibrateData(original, ammeter_scale=1e-6):
238 data = copy.deepcopy(original)
239 for i in range(0, len(data)):
240 data[i][1] = 16.8 * float(data[i][1]) / 4000.0
241 data[i][2] = ammeter_scale * 0.170 * float(data[i][2]) / 268.0
242 data[i][3] = ammeter_scale * 0.170 * float(data[i][3]) / 268.0
245 def ShowTCS(filename, raw=True,calibrate=True, normalise=False, show_error=False, plot=gnuplot.plot,with_="lp", step=1, output=None, title="", master_title="", smooth=0, show_peak=False, inflection=1):
251 if type(filename) == type(""):
252 data = GetData(filename)
255 filename = "tcs data"
258 title = BaseName(filename)
264 if type(smooth) == type([]):
265 for i in range(smooth[0]):
266 data = Smooth(data, m=smooth[1])
268 data = Smooth(data, m=smooth)
272 data = CalibrateData(data)
273 units = ["V", "uA / V"]
275 units = ["DAC counts", "ADC counts / DAC counts"]
278 gnuplot("set ylabel \"dI(E)/dE ("+str(units[1])+")\"")
280 data = MaxNormalise(data)
281 gnuplot("set ylabel \"dI(E)/dE (normalised)\"")
283 if (output != None and type(output) == type("")):
284 gnuplot("set term png size 640,480")
285 gnuplot("set output \""+str(output)+"\"")
287 if master_title == "":
288 master_title = "Total Current Spectrum S(E)"
289 if type(filename) == type("") and plot == gnuplot.plot:
290 if filename != "tcs data":
291 p = ReadParameters(filename)
293 master_title += "\\nSample: "+p["Sample"]
295 gnuplot("set title \""+str(master_title)+"\"")
296 gnuplot("set xlabel \"U ("+str(units[0])+")\"")
300 d = Derivative(data, 1, 2, step=step)
304 ymax = 0.01 + 1.2 * max(d, key=lambda e : e[2])[2]
305 ymin = -0.01 + 1.2 * min(d, key=lambda e : e[2])[2]
306 gnuplot("set yrange ["+str(ymin)+":"+str(ymax)+"]")
309 plotList.append(Gnuplot.Data(d, using="2:3", with_=with_,title=title))
312 error1 = Derivative(data, 1, 2, -3,step=step)
313 error2 = Derivative(data, 1, 2, +3,step=step)
314 plotList.append(Gnuplot.Data(error1, using="2:3", with_=w,title="-sigma/2"))
315 plotList.append(Gnuplot.Data(error2, using="2:3", with_=w, title="+sigma/2"))
318 peak = SmoothPeakFind(d, ap=DoNothing, stop=1, inflection=inflection)
319 plotList += PlotPeaks(peak,with_="l lt -1", plot=None)
327 if (output != None and type(output) == type("")):
328 gnuplot("set term wxt")
334 def ShowData(filename,calibrate=True, normalise=False, show_error=False, plot=gnuplot.plot,with_="lp", step=1, output=None, title="", master_title="Sample Current I(E)", smooth=0):
335 if type(filename) == type(""):
336 data = GetData(filename)
339 filename = "raw data"
342 title = BaseName(filename)
350 if type(data) == type([]):
351 for i in range(0, smooth[0]):
352 data = Smooth(data, m=smooth[1])
354 data = Smooth(data, m = smooth)
357 data = CalibrateData(data)
360 units = ["DAC counts", "ADC counts"]
363 gnuplot("set ylabel \"I(E) ("+str(units[1])+")\"")
365 data = MaxNormalise(data)
366 gnuplot("set ylabel \"I(E) (normalised)\"")
368 if (output != None and type(output) == type("")):
369 gnuplot("set term png size 640,480")
370 gnuplot("set output \""+str(output)+"\"")
372 gnuplot("set title \""+str(master_title)+"\"")
373 gnuplot("set xlabel \"U ("+str(units[0])+")\"")
376 ymax = 0.005 + 1.2 * max(d, key=lambda e : e[2])[2]
377 ymin = -0.005 + 1.2 * min(d, key=lambda e : e[2])[2]
378 gnuplot("set yrange ["+str(ymin)+":"+str(ymax)+"]")
380 #d = Derivative(data, 1, 2, step=step)
384 plotList.append(Gnuplot.Data(data, using="2:3", with_=with_,title=title))
387 error1 = copy.deepcopy(data)
388 error2 = copy.deepcopy(data)
389 for i in range(len(data)):
391 error1[i][2] -= 0.50*float(data[i][3])
392 error2[i][2] += 0.50*float(data[i][3])
393 plotList.append(Gnuplot.Data(error1, using="2:3", with_=w,title="Error : Low bound"))
394 plotList.append(Gnuplot.Data(error2, using="2:3", with_=w, title="Error : Upper bound"))
399 if (output != None and type(output) == type("")):
400 gnuplot("set term wxt")
405 def ReadParameters(filename):
406 parameters = odict.odict()
407 input_file = open(filename, "r")
408 for line in input_file:
413 item = k[0].strip("# \r\n")
414 value = k[1].strip("# \r\n")
415 if (item in parameters):
416 parameters[item] = value
418 parameters.update({str(item) : value})
422 def PlotParameters(filename):
423 ReadParameters(filename)
425 def Smooth(data, m, k=2):
426 smooth = copy.deepcopy(data)
427 for i in range(len(smooth)):
430 for j in range(i-m,i+m):
431 if j >= 0 and j < len(smooth):
433 smooth[i][k] += data[j][k]
435 smooth[i][k] = smooth[i][k] / float(count)
437 smooth[i][k] = data[i][k]
441 def PeakFind(data, k=2,threshold=0.00, inflection=0):
443 for i in range(len(data)):
444 if i == 0 or i == len(data)-1:
446 #if abs(data[i][k]) < threshold * abs(max(data, key = lambda e : abs(e[k]))[k]):
449 left = data[i-1][k] - data[i][k]
450 right = data[i+1][k] - data[i][k]
451 if abs(left) < threshold*abs(data[i][k]):
453 if abs(right) < threshold*abs(data[i][k]):
456 results.append(data[i] + [inflection])
459 results += PeakFind(Derivative(data), k=k, threshold=threshold, inflection=inflection-1)
463 def SmoothPeakFind(data, a=1, k=2, ap=DoNothing, stop=10,smooth=5, inflection=0):
472 peaks = PeakFind(ap(s),k=k, inflection=inflection)
473 #print "m = " +str(m)
476 [add.append(f) for f in p]
479 #print "*New peak at " + str(p)
480 peakList.append([add])
483 for i in range(len(peakList)):
484 p2 = peakList[i][len(peakList[i])-1]
487 score.append([i, abs(p[a] - p2[1+a])])
489 score.sort(key = lambda e : e[1])
490 if len(score) == 0 or score[0][1] > 100:
491 #print "New peak at " + str(p)
492 peakList.append([add])
494 #print "Peak exists near " + str(p) + " ("+str(score[0][1])+") " + str(peakList[score[0][0]][len(peakList[score[0][0]])-1])
495 peakList[score[0][0]].append(add)
499 #results.append([m, []])
500 #[results[len(results)-1].append(f) for f in p]
502 s = Smooth(s, m=smooth,k=k)
504 #results.sort(key = lambda e : e[2])
511 def PlotPeaks(peaks, calibrate=True, with_="lp", plot=gnuplot.replot):
515 p.append(copy.deepcopy(p[len(p)-1]))
520 #print "Adding " + str(p) + " to list"
524 with_ = with_.split(" lt")[0] + " lt 9"
525 plotList.append(Gnuplot.Data(p, using="3:1", with_=with_))
529 if len(plotList) > 0 and plot != None:
540 if (len(sys.argv) < 2):
541 sys.stderr.write(sys.argv[0] + " - Require arguments (filename)\n")
551 while i < len(sys.argv):
552 if sys.argv[i] == "--raw":
554 elif sys.argv[i] == "--tcs":
555 plotFunc = ShowTCS #lambda e : ShowTCS(e, show_peak=False)
556 elif sys.argv[i] == "--output":
557 if i+1 >= len(sys.argv):
558 sys.stderr.write("Need argument for "+sys.argv[i]+" switch\n")
560 gnuplot("set term postscript colour")
561 gnuplot("set output \""+sys.argv[i+1]+"\"")
563 elif sys.argv[i] == "--wxt":
564 gnuplot("set term wxt")
565 elif sys.argv[i] == "--normalise":
567 elif sys.argv[i] == "--unnormalise":
569 elif sys.argv[i] == "--title":
570 if i+1 >= len(sys.argv):
571 sys.stderr.write("Need argument for "+sys.argv[i]+" switch\n")
573 title = sys.argv[i+1]
575 elif sys.argv[i] == "--master_title":
576 if i+1 >= len(sys.argv):
577 sys.stderr.write("Need argument for "+sys.argv[i]+" switch\n")
579 master_title = sys.argv[i+1]
581 elif sys.argv[i] == "--smooth":
582 if i+1 >= len(sys.argv):
583 sys.stderr.write("Need argument for "+sys.argv[i]+" switch\n")
585 smooth = sys.argv[i+1]
586 smooth = map(int, smooth.split("x"))
590 elif sys.argv[i] == "--with":
591 if i+1 >= len(sys.argv):
592 sys.stderr.write("Need argument for "+sys.argv[i]+" switch\n")
594 with_ = sys.argv[i+1]
596 elif sys.argv[i] == "--output":
597 if i+1 >= len(sys.argv):
598 sys.stderr.write("Need argument for "+sys.argv[i]+" switch\n")
600 gnuplot("set term postscript colour")
601 gnuplot("set output \""+str(argv[i+1])+"\"")
603 plotFunc(sys.argv[i], plot=gnuplot.replot, normalise=normalise, title=title, master_title=master_title, smooth=smooth, with_=with_)
607 print "Done. Press enter to exit, or type name of file to save as."
608 out = sys.stdin.readline().strip("\t\r\n #")
610 gnuplot("set term postscript colour")
611 gnuplot("set output \""+out+"\"")
615 def ModelTCS(f, sigma, Emin, Emax, dE):
619 S = (1 - sigma(0))*f(-E) + FuncIntegrate(lambda e : f(e - E) * FuncDerivative(sigma, E, dE), Emin, Emax, dE)
620 data.append([0.00, E, S,0.00])
624 def IntegrateTCS(data, imin, imax=0, di=1):
631 total += data[i][2] * (data[i+1][1] - data[i][1])
635 def FuncIntegrate(f, xmin, xmax, dx):
643 def FuncDerivative(f, x, dx):
644 return 0.50*(f(x+dx) - f(x-dx))/dx
646 def FitTCS(data, min_mse=1e-4, max_fail=100, max_adjust=4,divide=10, plot=gnuplot.plot,smooth=0):
647 if type(data) == type(""):
656 if type(smooth) == type([]):
657 for _ in range(smooth[0]):
658 d = Smooth(d, m=smooth[1])
660 d = Smooth(d, m=smooth)
664 plotItems = ShowTCS(d, raw=False,smooth=smooth,plot=None)
665 plotItems.append(None)
667 peaks = SmoothPeakFind(d, smooth=5, stop=1, inflection=0)
668 peaks.sort(key = lambda e : e[len(e)-1][1])
672 for i in range(0,len(peaks)):
677 fits.append([l[3], l[2], 1.0])
680 l = peaks[i-2][len(peaks[i-2])]
682 fits.append([l[3], l[2], 1.0])
683 if i+2 <= len(peaks)-1:
684 l = peaks[i+2][len(peaks[i+2])]
686 fits.append([l[3], l[2], 1.0])
688 for i in range(len(fits)):
692 left = fits[i-1][1] - fits[i][1]
694 right = fits[i+1][1] - fits[i][1]
696 fits[i][2] = min([abs(0.5*left), abs(0.5*right)])
699 #print "Fits are " + str(fits)
707 dt = f[0] * gaussian(E - f[1], f[2])
708 #print " Increase total by " + str(dt)
710 #print "tcs returns " + str(total)
723 while failcount < max_fail and mse > min_mse:
724 i = random.randint(0, len(fits)-1)
725 j = random.randint(0, len(fits[i])-1)
727 # j = random.randint(0, len(fits[i])-1)
729 #print "Adjust " + str(i) + ","+str(j) + ": Iteration " + str(iterations) + " mse: " + str(mse)
734 fits[i][j] += adjust * (random.random() - 0.50)
736 while fits[i][j] <= 0.0005:
737 fits[i][j] = adjust * (random.random() - 0.50)
740 model = table(lambda e : [0.00, e, tcs(e), 0.00], 0, 16.8, divide*d[len(d)-1][1]/(len(d)))
741 mse = MeanSquareError(model, d[0::divide])
745 if failcount > max_fail / 2:
746 if adjust > 1.0/(2.0**max_adjust):
757 #model = table(lambda e : [0.00, e, tcs(e), 0.00], 0, 16.8, 16.8/len(d))
758 plotItems[len(plotItems)-1] = Gnuplot.Data(model, using="2:3", with_="l lt 3", title="model")
761 fits.sort(key = lambda e : e[0] * gaussian(0, e[2]), reverse=True)
764 gnuplot("set title \"MSE = "+str(mse)+"\\nfailcount = "+str(failcount)+"\\nadjust = "+str(adjust)+"\"")
765 gnuplot.plot(*plotItems)
769 return [fits, model,plotItems]
774 def SaveFit(filename, fit):
775 out = open(filename, "w", 0)
776 out.write("# TCS Fit\n")
779 out.write(str(f[0]) + "\t" + str(f[1]) + "\t" + str(f[2]) + "\n")
784 def LoadFit(filename):
785 infile = open(filename, "r")
786 if (infile.readline() != "# TCS Fit\n"):
787 sys.stderr.write("Error loading fit from file " + str(filename) + "\n")
792 f = infile.readline().strip("# \r\n\t").split("\t")
795 fit.append(map(float, f))
798 fit.sort(key = lambda e : e[0] * gaussian(0, e[2]), reverse=True)
802 # total += f[0] * gaussian(e - f[1], f[2])
805 #return table(lambda e : [0.0, e, model(e), 0.0], 0.0, 16.8, 16.8/400)
808 def MeanSquareError(model, real, k = 2):
810 for i in range(len(real)):
812 mse += (model[i][k] - real[i][k])**2
823 def table(f, xmin, xmax, dx):
831 def gaussian(x, sigma):
834 return math.exp(- (x**2.0)/(2.0 * sigma**2.0)) / (sigma * (2.0 * math.pi)**0.50)
836 def step(x, sigma, T):
839 return 1.0 / (math.exp((x - sigma)/T) + 1.0)
841 if __name__ == "__main__":