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)
353 data = CalibrateData(data)
356 units = ["DAC counts", "ADC counts"]
359 gnuplot("set ylabel \"I(E) ("+str(units[1])+")\"")
361 data = MaxNormalise(data)
362 gnuplot("set ylabel \"I(E) (normalised)\"")
364 if (output != None and type(output) == type("")):
365 gnuplot("set term png size 640,480")
366 gnuplot("set output \""+str(output)+"\"")
368 gnuplot("set title \""+str(master_title)+"\"")
369 gnuplot("set xlabel \"U ("+str(units[0])+")\"")
372 #d = Derivative(data, 1, 2, step=step)
376 plotList.append(Gnuplot.Data(data, using="2:3", with_=with_,title=title))
379 error1 = copy.deepcopy(data)
380 error2 = copy.deepcopy(data)
381 for i in range(len(data)):
383 error1[i][2] -= 0.50*float(data[i][3])
384 error2[i][2] += 0.50*float(data[i][3])
385 plotList.append(Gnuplot.Data(error1, using="2:3", with_=w,title="Error : Low bound"))
386 plotList.append(Gnuplot.Data(error2, using="2:3", with_=w, title="Error : Upper bound"))
391 if (output != None and type(output) == type("")):
392 gnuplot("set term wxt")
397 def ReadParameters(filename):
398 parameters = odict.odict()
399 input_file = open(filename, "r")
400 for line in input_file:
405 item = k[0].strip("# \r\n")
406 value = k[1].strip("# \r\n")
407 if (item in parameters):
408 parameters[item] = value
410 parameters.update({str(item) : value})
414 def PlotParameters(filename):
415 ReadParameters(filename)
417 def Smooth(data, m, k=2):
418 smooth = copy.deepcopy(data)
419 for i in range(len(smooth)):
422 for j in range(i-m,i+m):
423 if j >= 0 and j < len(smooth):
425 smooth[i][k] += data[j][k]
427 smooth[i][k] = smooth[i][k] / float(count)
429 smooth[i][k] = data[i][k]
433 def PeakFind(data, k=2,threshold=0.00, inflection=0):
435 for i in range(len(data)):
436 if i == 0 or i == len(data)-1:
438 #if abs(data[i][k]) < threshold * abs(max(data, key = lambda e : abs(e[k]))[k]):
441 left = data[i-1][k] - data[i][k]
442 right = data[i+1][k] - data[i][k]
443 if abs(left) < threshold*abs(data[i][k]):
445 if abs(right) < threshold*abs(data[i][k]):
448 results.append(data[i] + [inflection])
451 results += PeakFind(Derivative(data), k=k, threshold=threshold, inflection=inflection-1)
455 def SmoothPeakFind(data, a=1, k=2, ap=DoNothing, stop=10,smooth=5, inflection=0):
464 peaks = PeakFind(ap(s),k=k, inflection=inflection)
465 #print "m = " +str(m)
468 [add.append(f) for f in p]
471 #print "*New peak at " + str(p)
472 peakList.append([add])
475 for i in range(len(peakList)):
476 p2 = peakList[i][len(peakList[i])-1]
479 score.append([i, abs(p[a] - p2[1+a])])
481 score.sort(key = lambda e : e[1])
482 if len(score) == 0 or score[0][1] > 100:
483 #print "New peak at " + str(p)
484 peakList.append([add])
486 #print "Peak exists near " + str(p) + " ("+str(score[0][1])+") " + str(peakList[score[0][0]][len(peakList[score[0][0]])-1])
487 peakList[score[0][0]].append(add)
491 #results.append([m, []])
492 #[results[len(results)-1].append(f) for f in p]
494 s = Smooth(s, m=smooth,k=k)
496 #results.sort(key = lambda e : e[2])
503 def PlotPeaks(peaks, calibrate=True, with_="lp", plot=gnuplot.replot):
507 p.append(copy.deepcopy(p[len(p)-1]))
512 #print "Adding " + str(p) + " to list"
516 with_ = with_.split(" lt")[0] + " lt 9"
517 plotList.append(Gnuplot.Data(p, using="3:1", with_=with_))
521 if len(plotList) > 0 and plot != None:
532 if (len(sys.argv) < 2):
533 sys.stderr.write(sys.argv[0] + " - Require arguments (filename)\n")
541 while i < len(sys.argv):
542 if sys.argv[i] == "--raw":
544 elif sys.argv[i] == "--tcs":
545 plotFunc = ShowTCS #lambda e : ShowTCS(e, show_peak=False)
546 elif sys.argv[i] == "--output":
547 if i+1 >= len(sys.argv):
548 sys.stderr.write("Need argument for "+sys.argv[i]+" switch\n")
550 gnuplot("set term postscript colour")
551 gnuplot("set output \""+sys.argv[i+1]+"\"")
553 elif sys.argv[i] == "--wxt":
554 gnuplot("set term wxt")
555 elif sys.argv[i] == "--normalise":
557 elif sys.argv[i] == "--unnormalise":
559 elif sys.argv[i] == "--title":
560 if i+1 >= len(sys.argv):
561 sys.stderr.write("Need argument for "+sys.argv[i]+" switch\n")
563 title = sys.argv[i+1]
565 elif sys.argv[i] == "--master_title":
566 if i+1 >= len(sys.argv):
567 sys.stderr.write("Need argument for "+sys.argv[i]+" switch\n")
569 master_title = sys.argv[i+1]
572 plotFunc(sys.argv[i], plot=gnuplot.replot, normalise=normalise, title=title, master_title=master_title)
576 print "Done. Press enter to exit, or type name of file to save as."
577 out = sys.stdin.readline().strip("\t\r\n #")
579 gnuplot("set term postscript colour")
580 gnuplot("set output \""+out+"\"")
584 def ModelTCS(f, sigma, Emin, Emax, dE):
588 S = (1 - sigma(0))*f(-E) + FuncIntegrate(lambda e : f(e - E) * FuncDerivative(sigma, E, dE), Emin, Emax, dE)
589 data.append([0.00, E, S,0.00])
593 def IntegrateTCS(data, imin, imax=0, di=1):
600 total += data[i][2] * (data[i+1][1] - data[i][1])
604 def FuncIntegrate(f, xmin, xmax, dx):
612 def FuncDerivative(f, x, dx):
613 return 0.50*(f(x+dx) - f(x-dx))/dx
615 def FitTCS(data, min_mse=1e-4, max_fail=100, max_adjust=4,divide=10, plot=gnuplot.plot,smooth=0):
616 if type(data) == type(""):
625 if type(smooth) == type([]):
626 for _ in range(smooth[0]):
627 d = Smooth(d, m=smooth[1])
629 d = Smooth(d, m=smooth)
633 plotItems = ShowTCS(d, raw=False,smooth=smooth,plot=None)
634 plotItems.append(None)
636 peaks = SmoothPeakFind(d, smooth=5, stop=1, inflection=0)
637 peaks.sort(key = lambda e : e[len(e)-1][1])
641 for i in range(0,len(peaks)):
646 fits.append([l[3], l[2], 1.0])
649 l = peaks[i-2][len(peaks[i-2])]
651 fits.append([l[3], l[2], 1.0])
652 if i+2 <= len(peaks)-1:
653 l = peaks[i+2][len(peaks[i+2])]
655 fits.append([l[3], l[2], 1.0])
657 for i in range(len(fits)):
661 left = fits[i-1][1] - fits[i][1]
663 right = fits[i+1][1] - fits[i][1]
665 fits[i][2] = min([abs(0.5*left), abs(0.5*right)])
668 #print "Fits are " + str(fits)
676 dt = f[0] * gaussian(E - f[1], f[2])
677 #print " Increase total by " + str(dt)
679 #print "tcs returns " + str(total)
692 while failcount < max_fail and mse > min_mse:
693 i = random.randint(0, len(fits)-1)
694 j = random.randint(0, len(fits[i])-1)
696 # j = random.randint(0, len(fits[i])-1)
698 #print "Adjust " + str(i) + ","+str(j) + ": Iteration " + str(iterations) + " mse: " + str(mse)
703 fits[i][j] += adjust * (random.random() - 0.50)
705 while fits[i][j] <= 0.0005:
706 fits[i][j] = adjust * (random.random() - 0.50)
709 model = table(lambda e : [0.00, e, tcs(e), 0.00], 0, 16.8, divide*d[len(d)-1][1]/(len(d)))
710 mse = MeanSquareError(model, d[0::divide])
714 if failcount > max_fail / 2:
715 if adjust > 1.0/(2.0**max_adjust):
726 #model = table(lambda e : [0.00, e, tcs(e), 0.00], 0, 16.8, 16.8/len(d))
727 plotItems[len(plotItems)-1] = Gnuplot.Data(model, using="2:3", with_="l lt 3", title="model")
730 fits.sort(key = lambda e : e[0] * gaussian(0, e[2]), reverse=True)
733 gnuplot("set title \"MSE = "+str(mse)+"\\nfailcount = "+str(failcount)+"\\nadjust = "+str(adjust)+"\"")
734 gnuplot.plot(*plotItems)
738 return [fits, model,plotItems]
743 def SaveFit(filename, fit):
744 out = open(filename, "w", 0)
745 out.write("# TCS Fit\n")
748 out.write(str(f[0]) + "\t" + str(f[1]) + "\t" + str(f[2]) + "\n")
753 def LoadFit(filename):
754 infile = open(filename, "r")
755 if (infile.readline() != "# TCS Fit\n"):
756 sys.stderr.write("Error loading fit from file " + str(filename) + "\n")
761 f = infile.readline().strip("# \r\n\t").split("\t")
764 fit.append(map(float, f))
767 fit.sort(key = lambda e : e[0] * gaussian(0, e[2]), reverse=True)
771 # total += f[0] * gaussian(e - f[1], f[2])
774 #return table(lambda e : [0.0, e, model(e), 0.0], 0.0, 16.8, 16.8/400)
777 def MeanSquareError(model, real, k = 2):
779 for i in range(len(real)):
781 mse += (model[i][k] - real[i][k])**2
792 def table(f, xmin, xmax, dx):
800 def gaussian(x, sigma):
803 return math.exp(- (x**2.0)/(2.0 * sigma**2.0)) / (sigma * (2.0 * math.pi)**0.50)
805 def step(x, sigma, T):
808 return 1.0 / (math.exp((x - sigma)/T) + 1.0)
810 if __name__ == "__main__":