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):
108 #print "Derivative called"
114 for i in range(0, len(data),step):
115 result[len(result)-1] = [d for d in data[i]]
117 dE[0] = data[i][a] - data[i-step][a]
118 dI[0] = data[i][b] - data[i-step][b]
122 if (i < len(data)-step):
123 dE[1] = data[i+step][a] - data[i][a]
124 dI[1] = data[i+step][b] - data[i][b]
129 #print str(data[i]) + " ["+str(sigma)+"] = " + str(data[i][int(abs(sigma))])
131 if dI[0] != None: dI[0] -= 0.5*data[i][int(abs(sigma))]
132 if dI[1] != None: dI[1] -= 0.5*data[i][int(abs(sigma))]
134 if dI[0] != None: dI[0] += 0.5*data[i][int(abs(sigma))]
135 if dI[1] != None: dI[1] += 0.5*data[i][int(abs(sigma))]
150 deltaI /= float(count)
151 deltaE /= float(count)
155 result[len(result)-1][b] = (deltaI / deltaE)
157 result[len(result)-1][b] = 0.0
159 result[len(result)-1][b] = 0.0
162 return result[0:len(result)-1]
164 def MaxNormalise(data, u=2):
165 result = copy.deepcopy(data)
168 maxval = max(data, key = lambda e : e[u])[u]
178 def Average(data_sets, u=1):
183 #print "Already have " + str(p[u])
185 for i in range(0, len(p)):
186 avg[p[u]][0][i] += p[i]
188 #print "Create key for " + str(p[u])
189 avg.update({p[u] : [p, 1]})
192 for i in range(0, len(avg[a][0])):
193 avg[a][0][i] /= float(avg[a][1])
194 return map(lambda e : e[1][0], sorted(avg.items(), key = lambda e : e[0]))
196 def FullWidthAtHalfMax(data, u=1):
197 maxval = max(data, key = lambda e : e[u])
198 peak = data.index(maxval)
202 for i in range(1, len(data)/2):
204 if (peak-i > 0 and data[peak-i] < 0.50*maxval):
205 lhs = data[peak-i][u]
207 if (peak+i < len(data) and data[peak+i] < 0.50*maxval):
209 if lhs != None and rhs != None:
211 if rhs == None or lhs == None:
212 return abs(data[len(data)-1][0] - data[0][0])
214 return abs(rhs - lhs)
216 def SaveData(filename, data):
217 out = open(filename, "w", 0)
219 for i in range(0, len(a)):
225 def AverageAllData(directory, save=None, normalise=True):
227 if save == None: save = directory+"/average.dat"
228 for f in FindDataFiles(directory):
234 a = Average(data_sets)
238 def CalibrateData(original, ammeter_scale=1e-6):
239 data = copy.deepcopy(original)
240 for i in range(0, len(data)):
241 data[i][1] = 16.8 * float(data[i][1]) / 4000.0
242 data[i][2] = ammeter_scale * 0.170 * float(data[i][2]) / 268.0
243 data[i][3] = ammeter_scale * 0.170 * float(data[i][3]) / 268.0
246 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):
252 if type(filename) == type(""):
253 data = GetData(filename)
256 filename = "tcs data"
259 title = BaseName(filename)
265 if type(smooth) == type([]):
266 for i in range(smooth[0]):
267 data = Smooth(data, m=smooth[1])
269 data = Smooth(data, m=smooth)
273 data = CalibrateData(data)
274 units = ["Volts", "uA / V"]
276 units = ["DAC counts", "ADC counts / DAC counts"]
279 gnuplot("set ylabel \"S(U) ("+str(units[1])+")\"")
281 data = MaxNormalise(data)
282 gnuplot("set ylabel \"S(U) (arb. units)\"")
284 if (output != None and type(output) == type("")):
285 gnuplot("set term png size 640,480")
286 gnuplot("set output \""+str(output)+"\"")
288 if master_title == "":
289 master_title = "Total Current Spectrum S(U)"
290 if type(filename) == type("") and plot == gnuplot.plot:
291 if filename != "tcs data":
292 p = ReadParameters(filename)
294 master_title += "\\nSample: "+p["Sample"]
296 gnuplot("set title \""+str(master_title)+"\"")
297 gnuplot("set xlabel \"U ("+str(units[0])+")\"")
298 gnuplot("set xrange [0:16]")
299 gnuplot("set mxtics 10")
300 gnuplot("set mytics 10")
303 d = Derivative(data, 1, 2, step=step)
307 ymax = 0.01 + 1.2 * max(d, key=lambda e : e[2])[2]
308 ymin = -0.01 + 1.2 * min(d, key=lambda e : e[2])[2]
309 gnuplot("set yrange ["+str(ymin)+":"+str(ymax)+"]")
312 plotList.append(Gnuplot.Data(d, using="2:3", with_=with_,title=title))
315 error1 = Derivative(data, 1, 2, -3,step=step)
316 error2 = Derivative(data, 1, 2, +3,step=step)
317 plotList.append(Gnuplot.Data(error1, using="2:3", with_=w,title="-sigma/2"))
318 plotList.append(Gnuplot.Data(error2, using="2:3", with_=w, title="+sigma/2"))
321 peak = SmoothPeakFind(d, ap=DoNothing, stop=1, inflection=inflection)
322 plotList += PlotPeaks(peak,with_="l lt -1", plot=None)
330 if (output != None and type(output) == type("")):
331 gnuplot("set term wxt")
337 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):
338 if type(filename) == type(""):
339 data = GetData(filename)
342 filename = "raw data"
345 title = BaseName(filename)
353 if type(data) == type([]):
354 for i in range(0, smooth[0]):
355 data = Smooth(data, m=smooth[1])
357 data = Smooth(data, m = smooth)
360 data = CalibrateData(data)
363 units = ["DAC counts", "ADC counts"]
366 gnuplot("set ylabel \"I(U) ("+str(units[1])+")\"")
368 data = MaxNormalise(data)
369 gnuplot("set ylabel \"I(U) (arb. units)\"")
371 if (output != None and type(output) == type("")):
372 gnuplot("set term png size 640,480")
373 gnuplot("set output \""+str(output)+"\"")
375 gnuplot("set title \""+str(master_title)+"\"")
376 gnuplot("set xlabel \"U ("+str(units[0])+")\"")
379 ymax = 0.005 + 1.2 * max(data, key=lambda e : e[2])[2]
380 ymin = -0.005 + 1.2 * min(data, key=lambda e : e[2])[2]
381 gnuplot("set yrange ["+str(ymin)+":"+str(ymax)+"]")
383 #d = Derivative(data, 1, 2, step=step)
387 plotList.append(Gnuplot.Data(data, using="2:3", with_=with_,title=title))
390 error1 = copy.deepcopy(data)
391 error2 = copy.deepcopy(data)
392 for i in range(len(data)):
394 error1[i][2] -= 0.50*float(data[i][3])
395 error2[i][2] += 0.50*float(data[i][3])
396 plotList.append(Gnuplot.Data(error1, using="2:3", with_=w,title="Error : Low bound"))
397 plotList.append(Gnuplot.Data(error2, using="2:3", with_=w, title="Error : Upper bound"))
402 if (output != None and type(output) == type("")):
403 gnuplot("set term wxt")
408 def ReadParameters(filename):
409 parameters = odict.odict()
410 input_file = open(filename, "r")
411 for line in input_file:
416 item = k[0].strip("# \r\n")
417 value = k[1].strip("# \r\n")
418 if (item in parameters):
419 parameters[item] = value
421 parameters.update({str(item) : value})
425 def PlotParameters(filename):
426 ReadParameters(filename)
428 def Smooth(data, m, k=2):
429 smooth = copy.deepcopy(data)
430 for i in range(len(smooth)):
433 for j in range(i-m,i+m):
434 if j >= 0 and j < len(smooth):
436 smooth[i][k] += data[j][k]
438 smooth[i][k] = smooth[i][k] / float(count)
440 smooth[i][k] = data[i][k]
444 def PeakFind(data, k=2,threshold=0.00, inflection=0):
446 for i in range(len(data)):
447 if i == 0 or i == len(data)-1:
449 #if abs(data[i][k]) < threshold * abs(max(data, key = lambda e : abs(e[k]))[k]):
452 left = data[i-1][k] - data[i][k]
453 right = data[i+1][k] - data[i][k]
454 if abs(left) < threshold*abs(data[i][k]):
456 if abs(right) < threshold*abs(data[i][k]):
459 results.append(data[i] + [inflection])
462 results += PeakFind(Derivative(data), k=k, threshold=threshold, inflection=inflection-1)
466 def SmoothPeakFind(data, a=1, k=2, ap=DoNothing, stop=10,smooth=5, inflection=0):
475 peaks = PeakFind(ap(s),k=k, inflection=inflection)
476 #print "m = " +str(m)
479 [add.append(f) for f in p]
482 #print "*New peak at " + str(p)
483 peakList.append([add])
486 for i in range(len(peakList)):
487 p2 = peakList[i][len(peakList[i])-1]
490 score.append([i, abs(p[a] - p2[1+a])])
492 score.sort(key = lambda e : e[1])
493 if len(score) == 0 or score[0][1] > 100:
494 #print "New peak at " + str(p)
495 peakList.append([add])
497 #print "Peak exists near " + str(p) + " ("+str(score[0][1])+") " + str(peakList[score[0][0]][len(peakList[score[0][0]])-1])
498 peakList[score[0][0]].append(add)
502 #results.append([m, []])
503 #[results[len(results)-1].append(f) for f in p]
505 s = Smooth(s, m=smooth,k=k)
507 #results.sort(key = lambda e : e[2])
514 def PlotPeaks(peaks, calibrate=True, with_="lp", plot=gnuplot.replot):
518 p.append(copy.deepcopy(p[len(p)-1]))
523 #print "Adding " + str(p) + " to list"
527 with_ = with_.split(" lt")[0] + " lt 9"
528 plotList.append(Gnuplot.Data(p, using="3:1", with_=with_))
532 if len(plotList) > 0 and plot != None:
543 if (len(sys.argv) < 2):
544 sys.stderr.write(sys.argv[0] + " - Require arguments (filename)\n")
554 while i < len(sys.argv):
555 if sys.argv[i] == "--raw":
557 elif sys.argv[i] == "--tcs":
558 plotFunc = ShowTCS #lambda e : ShowTCS(e, show_peak=False)
559 elif sys.argv[i] == "--output":
560 if i+1 >= len(sys.argv):
561 sys.stderr.write("Need argument for "+sys.argv[i]+" switch\n")
563 gnuplot("set term postscript colour")
564 gnuplot("set output \""+sys.argv[i+1]+"\"")
566 elif sys.argv[i] == "--wxt":
567 gnuplot("set term wxt")
568 elif sys.argv[i] == "--normalise":
570 elif sys.argv[i] == "--unnormalise":
572 elif sys.argv[i] == "--title":
573 if i+1 >= len(sys.argv):
574 sys.stderr.write("Need argument for "+sys.argv[i]+" switch\n")
576 title = sys.argv[i+1]
578 elif sys.argv[i] == "--master_title":
579 if i+1 >= len(sys.argv):
580 sys.stderr.write("Need argument for "+sys.argv[i]+" switch\n")
582 master_title = sys.argv[i+1]
584 elif sys.argv[i] == "--smooth":
585 if i+1 >= len(sys.argv):
586 sys.stderr.write("Need argument for "+sys.argv[i]+" switch\n")
588 smooth = sys.argv[i+1]
589 smooth = map(int, smooth.split("x"))
593 elif sys.argv[i] == "--with":
594 if i+1 >= len(sys.argv):
595 sys.stderr.write("Need argument for "+sys.argv[i]+" switch\n")
597 with_ = sys.argv[i+1]
599 elif sys.argv[i] == "--output":
600 if i+1 >= len(sys.argv):
601 sys.stderr.write("Need argument for "+sys.argv[i]+" switch\n")
603 gnuplot("set term postscript colour")
604 gnuplot("set output \""+str(argv[i+1])+"\"")
606 plotFunc(sys.argv[i], plot=gnuplot.replot, normalise=normalise, title=title, master_title=master_title, smooth=smooth, with_=with_)
610 print "Done. Press enter to exit, or type name of file to save as."
611 out = sys.stdin.readline().strip("\t\r\n #")
613 gnuplot("set term postscript colour enhanced \"Arial Bold,18\"")
614 gnuplot("set output \""+out+"\"")
618 def ModelTCS(f, sigma, Emin, Emax, dE):
622 S = (1 - sigma(0))*f(-E) + FuncIntegrate(lambda e : f(e - E) * FuncDerivative(sigma, E, dE), Emin, Emax, dE)
623 data.append([0.00, E, S,0.00])
627 def IntegrateTCS(data, imin, imax=0, di=1):
634 total += data[i][2] * (data[i+1][1] - data[i][1])
638 def FuncIntegrate(f, xmin, xmax, dx):
646 def FuncDerivative(f, x, dx):
647 return 0.50*(f(x+dx) - f(x-dx))/dx
649 def FitTCS(data, min_mse=1e-4, max_fail=100, max_adjust=4,divide=10, plot=gnuplot.plot,smooth=0):
650 if type(data) == type(""):
659 if type(smooth) == type([]):
660 for _ in range(smooth[0]):
661 d = Smooth(d, m=smooth[1])
663 d = Smooth(d, m=smooth)
667 plotItems = ShowTCS(d, raw=False,smooth=smooth,plot=None)
668 plotItems.append(None)
670 peaks = SmoothPeakFind(d, smooth=5, stop=1, inflection=0)
671 peaks.sort(key = lambda e : e[len(e)-1][1])
675 for i in range(0,len(peaks)):
680 fits.append([l[3], l[2], 1.0])
683 l = peaks[i-2][len(peaks[i-2])]
685 fits.append([l[3], l[2], 1.0])
686 if i+2 <= len(peaks)-1:
687 l = peaks[i+2][len(peaks[i+2])]
689 fits.append([l[3], l[2], 1.0])
691 for i in range(len(fits)):
695 left = fits[i-1][1] - fits[i][1]
697 right = fits[i+1][1] - fits[i][1]
699 fits[i][2] = min([abs(0.5*left), abs(0.5*right)])
702 #print "Fits are " + str(fits)
710 dt = f[0] * gaussian(E - f[1], f[2])
711 #print " Increase total by " + str(dt)
713 #print "tcs returns " + str(total)
726 while failcount < max_fail and mse > min_mse:
727 i = random.randint(0, len(fits)-1)
728 j = random.randint(0, len(fits[i])-1)
730 # j = random.randint(0, len(fits[i])-1)
732 #print "Adjust " + str(i) + ","+str(j) + ": Iteration " + str(iterations) + " mse: " + str(mse)
737 fits[i][j] += adjust * (random.random() - 0.50)
739 while fits[i][j] <= 0.0005:
740 fits[i][j] = adjust * (random.random() - 0.50)
743 model = table(lambda e : [0.00, e, tcs(e), 0.00], 0, 16.8, divide*d[len(d)-1][1]/(len(d)))
744 mse = MeanSquareError(model, d[0::divide])
748 if failcount > max_fail / 2:
749 if adjust > 1.0/(2.0**max_adjust):
760 #model = table(lambda e : [0.00, e, tcs(e), 0.00], 0, 16.8, 16.8/len(d))
761 plotItems[len(plotItems)-1] = Gnuplot.Data(model, using="2:3", with_="l lt 3", title="model")
764 fits.sort(key = lambda e : e[0] * gaussian(0, e[2]), reverse=True)
767 gnuplot("set title \"MSE = "+str(mse)+"\\nfailcount = "+str(failcount)+"\\nadjust = "+str(adjust)+"\"")
768 gnuplot.plot(*plotItems)
772 return [fits, model,plotItems]
777 def SaveFit(filename, fit):
778 out = open(filename, "w", 0)
779 out.write("# TCS Fit\n")
782 out.write(str(f[0]) + "\t" + str(f[1]) + "\t" + str(f[2]) + "\n")
787 def LoadFit(filename):
788 infile = open(filename, "r")
789 if (infile.readline() != "# TCS Fit\n"):
790 sys.stderr.write("Error loading fit from file " + str(filename) + "\n")
795 f = infile.readline().strip("# \r\n\t").split("\t")
798 fit.append(map(float, f))
801 fit.sort(key = lambda e : e[0] * gaussian(0, e[2]), reverse=True)
805 # total += f[0] * gaussian(e - f[1], f[2])
808 #return table(lambda e : [0.0, e, model(e), 0.0], 0.0, 16.8, 16.8/400)
811 def MeanSquareError(model, real, k = 2):
813 for i in range(len(real)):
815 mse += (model[i][k] - real[i][k])**2
826 def table(f, xmin, xmax, dx):
834 def gaussian(x, sigma):
837 return math.exp(- (x**2.0)/(2.0 * sigma**2.0)) / (sigma * (2.0 * math.pi)**0.50)
839 def step(x, sigma, T):
842 return 1.0 / (math.exp((x - sigma)/T) + 1.0)
844 if __name__ == "__main__":