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):
55 input_file = open(filename, "r")
57 for line in input_file:
58 line = re.sub("#.*", "", line).strip("\r\n ")
62 line = map(lambda e : float(e), line.split("\t"))
65 for i in range(len(line)):
66 data[line[key]][0][i] += line[i]
67 data[line[key]][1] += 1
69 data.update({line[key] : [line, 1]})
71 d = map(lambda e : map(lambda f : float(f) / float(e[1][1]), e[1][0]), data.items())
72 d.sort(key = lambda e : e[key])
80 def AverageAllDataSets(directory=".", function=DoNothing):
82 for f in os.listdir(directory):
83 if os.path.isdir(directory+"/"+str(f)) == True:
85 for datafile in os.listdir(directory+"/"+str(f)):
86 if datafile.split(".")[1] == "dat":
87 data_set.append(GetData(f))
89 avg = Average(data_set)
90 dirs.update({f : avg})
93 def GetDataSets(directory="."):
95 for f in os.listdir(directory):
96 if os.path.isdir(directory+"/"+str(f)) == False:
97 if (len(f.split(".")) > 1 and f.split(".")[1] == "dat"):
98 d = GetData(directory+"/"+str(f))
105 def Derivative(data, a=1, b=2, sigma=None,step=1):
111 for i in range(0, len(data),step):
112 result[len(result)-1] = [d for d in data[i]]
114 dE[0] = data[i][a] - data[i-step][a]
115 dI[0] = data[i][b] - data[i-step][b]
119 if (i < len(data)-step):
120 dE[1] = data[i+step][a] - data[i][a]
121 dI[1] = data[i+step][b] - data[i][b]
126 #print str(data[i]) + " ["+str(sigma)+"] = " + str(data[i][int(abs(sigma))])
128 if dI[0] != None: dI[0] -= 0.5*data[i][int(abs(sigma))]
129 if dI[1] != None: dI[1] -= 0.5*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))]
147 deltaI /= float(count)
148 deltaE /= float(count)
152 result[len(result)-1][b] = (deltaI / deltaE)
154 result[len(result)-1][b] = 0.0
156 result[len(result)-1][b] = 0.0
159 return result[0:len(result)-1]
161 def MaxNormalise(data, u=2):
162 result = copy.deepcopy(data)
165 maxval = max(data, key = lambda e : e[u])[u]
172 def Average(data_sets, u=1):
177 #print "Already have " + str(p[u])
179 for i in range(0, len(p)):
180 avg[p[u]][0][i] += p[i]
182 #print "Create key for " + str(p[u])
183 avg.update({p[u] : [p, 1]})
186 for i in range(0, len(avg[a][0])):
187 avg[a][0][i] /= float(avg[a][1])
188 return map(lambda e : e[1][0], sorted(avg.items(), key = lambda e : e[0]))
190 def FullWidthAtHalfMax(data, u=1):
191 maxval = max(data, key = lambda e : e[u])
192 peak = data.index(maxval)
196 for i in range(1, len(data)/2):
198 if (peak-i > 0 and data[peak-i] < 0.50*maxval):
199 lhs = data[peak-i][u]
201 if (peak+i < len(data) and data[peak+i] < 0.50*maxval):
203 if lhs != None and rhs != None:
205 if rhs == None or lhs == None:
206 return abs(data[len(data)-1][0] - data[0][0])
208 return abs(rhs - lhs)
210 def SaveData(filename, data):
211 out = open(filename, "w", 0)
213 for i in range(0, len(a)):
219 def AverageAllData(directory, save=None):
221 if save == None: save = directory+"/average.dat"
222 for d in FindDataFiles(directory):
223 data_sets.append(GetData(d))
225 a = Average(data_sets)
229 def CalibrateData(original, ammeter_scale=1e-6):
230 data = copy.deepcopy(original)
231 for i in range(0, len(data)):
232 data[i][1] = 16.8 * float(data[i][1]) / 4000.0
233 data[i][2] = ammeter_scale * 0.170 * float(data[i][2]) / 268.0
234 data[i][3] = ammeter_scale * 0.170 * float(data[i][3]) / 268.0
237 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=True, inflection=1):
243 if type(filename) == type(""):
244 data = GetData(filename)
247 filename = "tcs data"
250 title = BaseName(filename)
256 if type(smooth) == type([]):
257 for i in range(smooth[0]):
258 data = Smooth(data, m=smooth[1])
260 data = Smooth(data, m=smooth)
264 data = CalibrateData(data)
265 units = ["V", "uA / V"]
267 units = ["DAC counts", "ADC counts / DAC counts"]
270 gnuplot("set ylabel \"dI(E)/dE ("+str(units[1])+")\"")
272 data = MaxNormalise(data)
273 gnuplot("set ylabel \"dI(E)/dE (normalised)\"")
275 if (output != None and type(output) == type("")):
276 gnuplot("set term png size 640,480")
277 gnuplot("set output \""+str(output)+"\"")
279 if master_title == "":
280 master_title = "Total Current Spectrum S(E)"
281 if type(filename) == type("") and plot == gnuplot.plot:
282 if filename != "tcs data":
283 p = ReadParameters(filename)
285 master_title += "\\nSample: "+p["Sample"]
287 gnuplot("set title \""+str(master_title)+"\"")
288 gnuplot("set xlabel \"U ("+str(units[0])+")\"")
292 d = Derivative(data, 1, 2, step=step)
296 ymax = 0.01 + 1.2 * max(d, key=lambda e : e[2])[2]
297 ymin = -0.01 + 1.2 * min(d, key=lambda e : e[2])[2]
298 gnuplot("set yrange ["+str(ymin)+":"+str(ymax)+"]")
301 plotList.append(Gnuplot.Data(d, using="2:3", with_=with_,title=title))
304 error1 = Derivative(data, 1, 2, -3,step=step)
305 error2 = Derivative(data, 1, 2, +3,step=step)
306 plotList.append(Gnuplot.Data(error1, using="2:3", with_=w,title="-sigma/2"))
307 plotList.append(Gnuplot.Data(error2, using="2:3", with_=w, title="+sigma/2"))
310 peak = SmoothPeakFind(d, ap=DoNothing, stop=1, inflection=inflection)
311 plotList += PlotPeaks(peak,with_="l lt -1", plot=None)
319 if (output != None and type(output) == type("")):
320 gnuplot("set term wxt")
326 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):
327 if type(filename) == type(""):
328 data = GetData(filename)
331 filename = "raw data"
334 title = BaseName(filename)
345 data = CalibrateData(data)
348 units = ["DAC counts", "ADC counts"]
351 gnuplot("set ylabel \"I(E) ("+str(units[1])+")\"")
353 data = MaxNormalise(data)
354 gnuplot("set ylabel \"I(E) (normalised)\"")
356 if (output != None and type(output) == type("")):
357 gnuplot("set term png size 640,480")
358 gnuplot("set output \""+str(output)+"\"")
360 gnuplot("set title \""+str(master_title)+"\"")
361 gnuplot("set xlabel \"U ("+str(units[0])+")\"")
364 #d = Derivative(data, 1, 2, step=step)
368 plotList.append(Gnuplot.Data(data, using="2:3", with_=with_,title=title))
371 error1 = copy.deepcopy(data)
372 error2 = copy.deepcopy(data)
373 for i in range(len(data)):
375 error1[i][2] -= 0.50*float(data[i][3])
376 error2[i][2] += 0.50*float(data[i][3])
377 plotList.append(Gnuplot.Data(error1, using="2:3", with_=w,title="Error : Low bound"))
378 plotList.append(Gnuplot.Data(error2, using="2:3", with_=w, title="Error : Upper bound"))
383 if (output != None and type(output) == type("")):
384 gnuplot("set term wxt")
389 def ReadParameters(filename):
390 parameters = odict.odict()
391 input_file = open(filename, "r")
392 for line in input_file:
397 item = k[0].strip("# \r\n")
398 value = k[1].strip("# \r\n")
399 if (item in parameters):
400 parameters[item] = value
402 parameters.update({str(item) : value})
406 def PlotParameters(filename):
407 ReadParameters(filename)
409 def Smooth(data, m, k=2):
410 smooth = copy.deepcopy(data)
411 for i in range(len(smooth)):
414 for j in range(i-m,i+m):
415 if j >= 0 and j < len(smooth):
417 smooth[i][k] += data[j][k]
419 smooth[i][k] = smooth[i][k] / float(count)
421 smooth[i][k] = data[i][k]
425 def PeakFind(data, k=2,threshold=0.00, inflection=0):
427 for i in range(len(data)):
428 if i == 0 or i == len(data)-1:
430 #if abs(data[i][k]) < threshold * abs(max(data, key = lambda e : abs(e[k]))[k]):
433 left = data[i-1][k] - data[i][k]
434 right = data[i+1][k] - data[i][k]
435 if abs(left) < threshold*abs(data[i][k]):
437 if abs(right) < threshold*abs(data[i][k]):
440 results.append(data[i] + [inflection])
443 results += PeakFind(Derivative(data), k=k, threshold=threshold, inflection=inflection-1)
447 def SmoothPeakFind(data, a=1, k=2, ap=DoNothing, stop=10,smooth=5, inflection=0):
456 peaks = PeakFind(ap(s),k=k, inflection=inflection)
457 #print "m = " +str(m)
460 [add.append(f) for f in p]
463 #print "*New peak at " + str(p)
464 peakList.append([add])
467 for i in range(len(peakList)):
468 p2 = peakList[i][len(peakList[i])-1]
471 score.append([i, abs(p[a] - p2[1+a])])
473 score.sort(key = lambda e : e[1])
474 if len(score) == 0 or score[0][1] > 100:
475 #print "New peak at " + str(p)
476 peakList.append([add])
478 #print "Peak exists near " + str(p) + " ("+str(score[0][1])+") " + str(peakList[score[0][0]][len(peakList[score[0][0]])-1])
479 peakList[score[0][0]].append(add)
483 #results.append([m, []])
484 #[results[len(results)-1].append(f) for f in p]
486 s = Smooth(s, m=smooth,k=k)
488 #results.sort(key = lambda e : e[2])
495 def PlotPeaks(peaks, calibrate=True, with_="lp", plot=gnuplot.replot):
499 p.append(copy.deepcopy(p[len(p)-1]))
504 #print "Adding " + str(p) + " to list"
508 with_ = with_.split(" lt")[0] + " lt 9"
509 plotList.append(Gnuplot.Data(p, using="3:1", with_=with_))
513 if len(plotList) > 0 and plot != None:
524 if (len(sys.argv) < 2):
525 sys.stderr.write(sys.argv[0] + " - Require arguments (filename)\n")
533 while i < len(sys.argv):
534 if sys.argv[i] == "--raw":
536 elif sys.argv[i] == "--tcs":
538 elif sys.argv[i] == "--output":
539 if i+1 >= len(sys.argv):
540 sys.stderr.write("Need argument for "+sys.argv[i]+" switch\n")
542 gnuplot("set term postscript colour")
543 gnuplot("set output \""+sys.argv[i+1]+"\"")
545 elif sys.argv[i] == "--wxt":
546 gnuplot("set term wxt")
547 elif sys.argv[i] == "--normalise":
549 elif sys.argv[i] == "--unnormalise":
551 elif sys.argv[i] == "--title":
552 if i+1 >= len(sys.argv):
553 sys.stderr.write("Need argument for "+sys.argv[i]+" switch\n")
555 title = sys.argv[i+1]
557 elif sys.argv[i] == "--master_title":
558 if i+1 >= len(sys.argv):
559 sys.stderr.write("Need argument for "+sys.argv[i]+" switch\n")
561 master_title = sys.argv[i+1]
564 plotFunc(sys.argv[i], plot=gnuplot.replot, normalise=normalise, title=title, master_title=master_title)
568 print "Done. Press enter to exit, or type name of file to save as."
569 out = sys.stdin.readline().strip("\t\r\n #")
571 gnuplot("set term postscript colour")
572 gnuplot("set output \""+out+"\"")
576 def ModelTCS(f, sigma, Emin, Emax, dE):
580 S = (1 - sigma(0))*f(-E) + FuncIntegrate(lambda e : f(e - E) * FuncDerivative(sigma, E, dE), Emin, Emax, dE)
581 data.append([0.00, E, S,0.00])
585 def IntegrateTCS(data, imin, imax=0, di=1):
592 total += data[i][2] * (data[i+1][1] - data[i][1])
596 def FuncIntegrate(f, xmin, xmax, dx):
604 def FuncDerivative(f, x, dx):
605 return 0.50*(f(x+dx) - f(x-dx))/dx
607 def FitTCS(data, min_mse=1e-4, max_fail=100, max_adjust=4,divide=10, plot=gnuplot.plot,smooth=0):
608 if type(data) == type(""):
617 if type(smooth) == type([]):
618 for _ in range(smooth[0]):
619 d = Smooth(d, m=smooth[1])
621 d = Smooth(d, m=smooth)
625 plotItems = ShowTCS(d, raw=False,smooth=smooth,plot=None)
626 plotItems.append(None)
628 peaks = SmoothPeakFind(d, smooth=5, stop=1, inflection=0)
629 peaks.sort(key = lambda e : e[len(e)-1][1])
633 for i in range(0,len(peaks)):
638 fits.append([l[3], l[2], 1.0])
641 l = peaks[i-2][len(peaks[i-2])]
643 fits.append([l[3], l[2], 1.0])
644 if i+2 <= len(peaks)-1:
645 l = peaks[i+2][len(peaks[i+2])]
647 fits.append([l[3], l[2], 1.0])
649 for i in range(len(fits)):
653 left = fits[i-1][1] - fits[i][1]
655 right = fits[i+1][1] - fits[i][1]
657 fits[i][2] = min([abs(0.5*left), abs(0.5*right)])
660 #print "Fits are " + str(fits)
668 dt = f[0] * gaussian(E - f[1], f[2])
669 #print " Increase total by " + str(dt)
671 #print "tcs returns " + str(total)
684 while failcount < max_fail and mse > min_mse:
685 i = random.randint(0, len(fits)-1)
686 j = random.randint(0, len(fits[i])-1)
688 # j = random.randint(0, len(fits[i])-1)
690 #print "Adjust " + str(i) + ","+str(j) + ": Iteration " + str(iterations) + " mse: " + str(mse)
695 fits[i][j] += adjust * (random.random() - 0.50)
697 while fits[i][j] <= 0.0005:
698 fits[i][j] = adjust * (random.random() - 0.50)
701 model = table(lambda e : [0.00, e, tcs(e), 0.00], 0, 16.8, divide*d[len(d)-1][1]/(len(d)))
702 mse = MeanSquareError(model, d[0::divide])
706 if failcount > max_fail / 2:
707 if adjust > 1.0/(2.0**max_adjust):
718 #model = table(lambda e : [0.00, e, tcs(e), 0.00], 0, 16.8, 16.8/len(d))
719 plotItems[len(plotItems)-1] = Gnuplot.Data(model, using="2:3", with_="l lt 3", title="model")
722 fits.sort(key = lambda e : e[0] * gaussian(0, e[2]), reverse=True)
725 gnuplot("set title \"MSE = "+str(mse)+"\\nfailcount = "+str(failcount)+"\\nadjust = "+str(adjust)+"\"")
726 gnuplot.plot(*plotItems)
730 return [fits, model,plotItems]
735 def SaveFit(filename, fit):
736 out = open(filename, "w", 0)
737 out.write("# TCS Fit\n")
740 out.write(str(f[0]) + "\t" + str(f[1]) + "\t" + str(f[2]) + "\n")
745 def LoadFit(filename):
746 infile = open(filename, "r")
747 if (infile.readline() != "# TCS Fit\n"):
748 sys.stderr.write("Error loading fit from file " + str(filename) + "\n")
753 f = infile.readline().strip("# \r\n\t").split("\t")
756 fit.append(map(float, f))
759 fit.sort(key = lambda e : e[0] * gaussian(0, e[2]), reverse=True)
763 # total += f[0] * gaussian(e - f[1], f[2])
766 #return table(lambda e : [0.0, e, model(e), 0.0], 0.0, 16.8, 16.8/400)
769 def MeanSquareError(model, real, k = 2):
771 for i in range(len(real)):
773 mse += (model[i][k] - real[i][k])**2
784 def table(f, xmin, xmax, dx):
792 def gaussian(x, sigma):
795 return math.exp(- (x**2.0)/(2.0 * sigma**2.0)) / (sigma * (2.0 * math.pi)**0.50)
797 def step(x, sigma, T):
800 return 1.0 / (math.exp((x - sigma)/T) + 1.0)
802 if __name__ == "__main__":