5 # @purpose Process TCS data
13 import re # Regular expressions - for removing comments
14 import odict #ordered dictionary
17 import Gnuplot, Gnuplot.funcutils
20 gnuplot = Gnuplot.Gnuplot()
23 gnuplot = Gnuplot.Gnuplot()
25 def FindDataFiles(directory=".", depth=1, result=None):
29 for f in os.listdir(directory):
30 if os.path.isdir(directory+"/"+str(f)):
32 result += FindDataFiles(directory+"/"+str(f), depth-1, result)
35 if (len(s) == 2 and s[1] == "dat"):
36 result.append(directory+"/"+str(f))
45 def DirectoryName(f, start=0,back=1):
47 return string.join(a[start:(len(a)-back)], "/")
49 def GetData(filename):
50 input_file = open(filename, "r")
52 for line in input_file:
53 line = re.sub("#.*", "", line).strip("\r\n ")
56 data.append(map(lambda e : float(e), line.split("\t")))
62 def AverageAllDataSets(directory=".", function=DoNothing):
64 for f in os.listdir(directory):
65 if os.path.isdir(directory+"/"+str(f)) == True:
67 for datafile in os.listdir(directory+"/"+str(f)):
68 if datafile.split(".")[1] == "dat":
69 data_set.append(GetData(f))
71 avg = Average(data_set)
72 dirs.update({f : avg})
75 def GetDataSets(directory="."):
77 for f in os.listdir(directory):
78 if os.path.isdir(directory+"/"+str(f)) == False:
79 if (len(f.split(".")) > 1 and f.split(".")[1] == "dat"):
80 d = GetData(directory+"/"+str(f))
87 def Derivative(data, a=1, b=2, sigma=None,step=1):
93 for i in range(0, len(data),step):
94 result[len(result)-1] = [d for d in data[i]]
96 dE[0] = data[i][a] - data[i-step][a]
97 dI[0] = data[i][b] - data[i-step][b]
101 if (i < len(data)-step):
102 dE[1] = data[i+step][a] - data[i][a]
103 dI[1] = data[i+step][b] - data[i][b]
108 #print str(data[i]) + " ["+str(sigma)+"] = " + str(data[i][int(abs(sigma))])
110 if dI[0] != None: dI[0] -= 0.5*data[i][int(abs(sigma))]
111 if dI[1] != None: dI[1] -= 0.5*data[i][int(abs(sigma))]
113 if dI[0] != None: dI[0] += 0.5*data[i][int(abs(sigma))]
114 if dI[1] != None: dI[1] += 0.5*data[i][int(abs(sigma))]
128 deltaI /= float(count)
129 deltaE /= float(count)
133 result[len(result)-1][b] = (deltaI / deltaE)
135 result[len(result)-1][b] = 0.0
138 return result[0:len(result)-1]
140 def MaxNormalise(data, u=2):
141 result = copy.deepcopy(data)
144 maxval = max(data, key = lambda e : e[u])[u]
151 def Average(data_sets, u=1):
156 #print "Already have " + str(p[u])
158 for i in range(0, len(p)):
159 avg[p[u]][0][i] += p[i]
161 #print "Create key for " + str(p[u])
162 avg.update({p[u] : [p, 1]})
165 for i in range(0, len(avg[a][0])):
166 avg[a][0][i] /= float(avg[a][1])
167 return map(lambda e : e[1][0], sorted(avg.items(), key = lambda e : e[0]))
169 def FullWidthAtHalfMax(data, u=1):
170 maxval = max(data, key = lambda e : e[u])
171 peak = data.index(maxval)
175 for i in range(1, len(data)/2):
177 if (peak-i > 0 and data[peak-i] < 0.50*maxval):
178 lhs = data[peak-i][u]
180 if (peak+i < len(data) and data[peak+i] < 0.50*maxval):
182 if lhs != None and rhs != None:
184 if rhs == None or lhs == None:
185 return abs(data[len(data)-1][0] - data[0][0])
187 return abs(rhs - lhs)
189 def SaveData(filename, data):
190 out = open(filename, "w", 0)
192 for i in range(0, len(a)):
198 def AverageAllData(directory, save=None)
200 if save == None: save = directory+"/average.dat":
201 for d in FindDataFiles(directory):
202 data_sets.append(GetData(d))
204 a = Average(data_sets)
208 def CalibrateData(original, ammeter_scale=1e-6):
209 data = copy.deepcopy(original)
210 for i in range(0, len(data)):
211 data[i][1] = 16.8 * float(data[i][1]) / 4000.0
212 data[i][2] = ammeter_scale * 0.170 * float(data[i][2]) / 268.0
213 data[i][3] = ammeter_scale * 0.170 * float(data[i][3]) / 268.0
216 def ShowTCS(filename, calibrate=True, normalise=False, show_error=False, plot=gnuplot.plot,with_="lp", step=1, output=None, title=""):
217 if type(filename) == type(""):
218 data = GetData(filename)
221 filename = "tcs data"
224 title = BaseName(filename)
230 data = CalibrateData(data)
231 units = ["V", "uA / V"]
233 units = ["DAC counts", "ADC counts / DAC counts"]
236 gnuplot("set ylabel \"dI(E)/dE ("+str(units[1])+")\"")
238 data = MaxNormalise(data)
239 gnuplot("set ylabel \"dI(E)/dE (normalised)\"")
241 if (output != None and type(output) == type("")):
242 gnuplot("set term png size 640,480")
243 gnuplot("set output \""+str(output)+"\"")
245 gnuplot("set title \"Total Current Spectrum S(E)\"")
246 gnuplot("set xlabel \"U ("+str(units[0])+")\"")
249 d = Derivative(data, 1, 2, step=step)
251 plot(Gnuplot.Data(d, using="2:3", with_=with_,title=title))
253 error1 = Derivative(data, 1, 2, -3,step=step)
254 error2 = Derivative(data, 1, 2, +3,step=step)
255 gnuplot.replot(Gnuplot.Data(error1, using="2:3", with_=w,title="-sigma/2"))
256 gnuplot.replot(Gnuplot.Data(error2, using="2:3", with_=w, title="+sigma/2"))
258 if (output != None and type(output) == type("")):
259 gnuplot("set term wxt")
262 def ShowData(filename,calibrate=True, normalise=False, show_error=False, plot=gnuplot.plot,with_="lp", step=1, output=None, title=""):
263 if type(filename) == type(""):
264 data = GetData(filename)
267 filename = "raw data"
270 title = BaseName(filename)
275 data = CalibrateData(data)
278 units = ["DAC counts", "ADC counts"]
281 gnuplot("set ylabel \"I(E) ("+str(units[1])+")\"")
283 data = MaxNormalise(data)
284 gnuplot("set ylabel \"I(E) (normalised)\"")
286 if (output != None and type(output) == type("")):
287 gnuplot("set term png size 640,480")
288 gnuplot("set output \""+str(output)+"\"")
290 gnuplot("set title \"Sample Current I(E)\"")
291 gnuplot("set xlabel \"U ("+str(units[0])+")\"")
294 #d = Derivative(data, 1, 2, step=step)
296 plot(Gnuplot.Data(data, using="2:3", with_=with_,title=title))
298 error1 = copy.deepcopy(data)
299 error2 = copy.deepcopy(data)
300 for i in range(len(data)):
302 error1[i][2] -= 0.50*float(data[i][3])
303 error2[i][2] += 0.50*float(data[i][3])
304 gnuplot.replot(Gnuplot.Data(error1, using="2:3", with_=w,title="Error : Low bound"))
305 gnuplot.replot(Gnuplot.Data(error2, using="2:3", with_=w, title="Error : Upper bound"))
307 if (output != None and type(output) == type("")):
308 gnuplot("set term wxt")
312 if (len(sys.argv) < 2):
313 sys.stderr.write(sys.argv[0] + " - Require arguments (filename)\n")
317 gnuplot("set style data lp")
318 gnuplot("set key outside right")
319 #gnuplot("set title \"Au on Si (50min 3.5A 3-6 e-8mbar)\"")
320 #gnuplot("set xlabel \"E (DAC Counts)\"")
321 #gnuplot("set ylabel \"S(E) (ADC/DAC Counts)\"")
322 #gnuplot("set term postscript colour")
323 #gnuplot("set output \"test.eps\"")
324 for i in range(1, len(sys.argv)):
325 tcs.append(Derivative(GetData(sys.argv[i]), 1, 2))
326 #tcs.append(GetTCS(GetData(sys.argv[i])))
327 if (len(tcs[i-1]) > 0):
328 gnuplot.replot(Gnuplot.Data(tcs[i-1], title=sys.argv[i], with_="lp"))
330 # Now average the data
336 sys.stdout.write(str(a[0]) + "\t" + str(a[1]) + "\t" + str(a[1]) + "\n")
337 gnuplot.replot(Gnuplot.Data(avg, title="Average", with_="l lw 2"))
339 sys.stdout.write("Save averaged data as (blank for no save): ")
340 filename = sys.stdin.readline().strip(" \r\n\t")
342 SaveData(filename, avg)
347 if __name__ == "__main__":