+def AverageAllData(directory, save=None, normalise=True):
+ data_sets = []
+ if save == None: save = directory+"/average.dat"
+ for f in FindDataFiles(directory):
+ d = GetData(f)
+ if normalise:
+ d = MaxNormalise(d)
+ data_sets.append(d)
+
+ a = Average(data_sets)
+ SaveData(save, a)
+ return a
+
+def CalibrateData(original, ammeter_scale=1e-6):
+ data = copy.deepcopy(original)
+ for i in range(0, len(data)):
+ data[i][1] = 16.8 * float(data[i][1]) / 4000.0
+ data[i][2] = ammeter_scale * 0.170 * float(data[i][2]) / 268.0
+ data[i][3] = ammeter_scale * 0.170 * float(data[i][3]) / 268.0
+ return data
+
+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):
+
+ if raw == False:
+ calibrate = False
+ normalise = False
+
+ if type(filename) == type(""):
+ data = GetData(filename)
+ else:
+ data = filename
+ filename = "tcs data"
+
+ if (title == ""):
+ title = BaseName(filename)
+
+ if (len(data) <= 0):
+ return data
+
+ if (smooth > 0):
+ if type(smooth) == type([]):
+ for i in range(smooth[0]):
+ data = Smooth(data, m=smooth[1])
+ else:
+ data = Smooth(data, m=smooth)
+
+
+ if calibrate:
+ data = CalibrateData(data)
+ units = ["V", "uA / V"]
+ else:
+ units = ["DAC counts", "ADC counts / DAC counts"]
+
+ if not normalise:
+ gnuplot("set ylabel \"dI(E)/dE ("+str(units[1])+")\"")
+ else:
+ data = MaxNormalise(data)
+ gnuplot("set ylabel \"dI(E)/dE (normalised)\"")
+
+ if (output != None and type(output) == type("")):
+ gnuplot("set term png size 640,480")
+ gnuplot("set output \""+str(output)+"\"")
+
+ if master_title == "":
+ master_title = "Total Current Spectrum S(E)"
+ if type(filename) == type("") and plot == gnuplot.plot:
+ if filename != "tcs data":
+ p = ReadParameters(filename)
+ if "Sample" in p:
+ master_title += "\\nSample: "+p["Sample"]
+
+ gnuplot("set title \""+str(master_title)+"\"")
+ gnuplot("set xlabel \"U ("+str(units[0])+")\"")
+
+
+ if raw:
+ d = Derivative(data, 1, 2, step=step)
+ else:
+ d = data
+
+ ymax = 0.01 + 1.2 * max(d, key=lambda e : e[2])[2]
+ ymin = -0.01 + 1.2 * min(d, key=lambda e : e[2])[2]
+ gnuplot("set yrange ["+str(ymin)+":"+str(ymax)+"]")
+
+ plotList = []
+ plotList.append(Gnuplot.Data(d, using="2:3", with_=with_,title=title))
+
+ if (show_error):
+ error1 = Derivative(data, 1, 2, -3,step=step)
+ error2 = Derivative(data, 1, 2, +3,step=step)
+ plotList.append(Gnuplot.Data(error1, using="2:3", with_=w,title="-sigma/2"))
+ plotList.append(Gnuplot.Data(error2, using="2:3", with_=w, title="+sigma/2"))
+
+ if (show_peak):
+ peak = SmoothPeakFind(d, ap=DoNothing, stop=1, inflection=inflection)
+ plotList += PlotPeaks(peak,with_="l lt -1", plot=None)
+
+
+
+ if (plot != None):
+ plot(*plotList)
+ time.sleep(0.2)
+
+ if (output != None and type(output) == type("")):
+ gnuplot("set term wxt")
+
+ if (plot == None):
+ return plotList
+ return data
+
+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):
+ if type(filename) == type(""):
+ data = GetData(filename)
+ else:
+ data = filename
+ filename = "raw data"
+
+ if (title == ""):
+ title = BaseName(filename)
+
+
+ if len(data) <= 0:
+ return data
+
+
+ if (smooth > 0):
+ data = Smooth(data)
+
+ if calibrate:
+ data = CalibrateData(data)
+ units = ["V", "uA"]
+ else:
+ units = ["DAC counts", "ADC counts"]
+
+ if not normalise:
+ gnuplot("set ylabel \"I(E) ("+str(units[1])+")\"")
+ else:
+ data = MaxNormalise(data)
+ gnuplot("set ylabel \"I(E) (normalised)\"")
+
+ if (output != None and type(output) == type("")):
+ gnuplot("set term png size 640,480")
+ gnuplot("set output \""+str(output)+"\"")
+
+ gnuplot("set title \""+str(master_title)+"\"")
+ gnuplot("set xlabel \"U ("+str(units[0])+")\"")
+
+
+ #d = Derivative(data, 1, 2, step=step)
+
+ plotList = []
+
+ plotList.append(Gnuplot.Data(data, using="2:3", with_=with_,title=title))
+ time.sleep(0.1)
+ if (show_error):
+ error1 = copy.deepcopy(data)
+ error2 = copy.deepcopy(data)
+ for i in range(len(data)):
+ #print str(data[i])
+ error1[i][2] -= 0.50*float(data[i][3])
+ error2[i][2] += 0.50*float(data[i][3])
+ plotList.append(Gnuplot.Data(error1, using="2:3", with_=w,title="Error : Low bound"))
+ plotList.append(Gnuplot.Data(error2, using="2:3", with_=w, title="Error : Upper bound"))
+
+ if plot != None:
+
+ plot(*plotList)
+ if (output != None and type(output) == type("")):
+ gnuplot("set term wxt")
+ return data
+ else:
+ return plotList
+
+def ReadParameters(filename):
+ parameters = odict.odict()
+ input_file = open(filename, "r")
+ for line in input_file:
+ k = line.split("=")
+ item = None
+ #print str(k)
+ if (len(k) >= 2):
+ item = k[0].strip("# \r\n")
+ value = k[1].strip("# \r\n")
+ if (item in parameters):
+ parameters[item] = value
+ else:
+ parameters.update({str(item) : value})
+ input_file.close()
+ return parameters
+
+def PlotParameters(filename):
+ ReadParameters(filename)
+
+def Smooth(data, m, k=2):
+ smooth = copy.deepcopy(data)
+ for i in range(len(smooth)):
+ count = 0
+ smooth[i][k] = 0.0
+ for j in range(i-m,i+m):
+ if j >= 0 and j < len(smooth):
+ count += 1
+ smooth[i][k] += data[j][k]
+ if count > 0:
+ smooth[i][k] = smooth[i][k] / float(count)
+ else:
+ smooth[i][k] = data[i][k]
+
+ return smooth
+
+def PeakFind(data, k=2,threshold=0.00, inflection=0):
+ results = []
+ for i in range(len(data)):
+ if i == 0 or i == len(data)-1:
+ continue
+ #if abs(data[i][k]) < threshold * abs(max(data, key = lambda e : abs(e[k]))[k]):
+ # continue
+
+ left = data[i-1][k] - data[i][k]
+ right = data[i+1][k] - data[i][k]
+ if abs(left) < threshold*abs(data[i][k]):
+ continue
+ if abs(right) < threshold*abs(data[i][k]):
+ continue
+ if left*right > 0:
+ results.append(data[i] + [inflection])
+
+ if inflection > 0:
+ results += PeakFind(Derivative(data), k=k, threshold=threshold, inflection=inflection-1)
+
+ return results
+
+def SmoothPeakFind(data, a=1, k=2, ap=DoNothing, stop=10,smooth=5, inflection=0):
+ s = data
+ #results = []
+
+ peakList = []
+
+ m = 0
+ while m < stop:
+ #results.append([])
+ peaks = PeakFind(ap(s),k=k, inflection=inflection)
+ #print "m = " +str(m)
+ for p in peaks:
+ add = [m]
+ [add.append(f) for f in p]
+
+ if m == 0:
+ #print "*New peak at " + str(p)
+ peakList.append([add])
+ else:
+ score = []
+ for i in range(len(peakList)):
+ p2 = peakList[i][len(peakList[i])-1]
+ if m - p2[0] > 1:
+ continue
+ score.append([i, abs(p[a] - p2[1+a])])
+
+ score.sort(key = lambda e : e[1])
+ if len(score) == 0 or score[0][1] > 100:
+ #print "New peak at " + str(p)
+ peakList.append([add])
+ else:
+ #print "Peak exists near " + str(p) + " ("+str(score[0][1])+") " + str(peakList[score[0][0]][len(peakList[score[0][0]])-1])
+ peakList[score[0][0]].append(add)
+
+
+
+ #results.append([m, []])
+ #[results[len(results)-1].append(f) for f in p]
+ m += 1
+ s = Smooth(s, m=smooth,k=k)
+
+ #results.sort(key = lambda e : e[2])
+
+ #peaks = []
+ return peakList
+
+
+
+def PlotPeaks(peaks, calibrate=True, with_="lp", plot=gnuplot.replot):
+
+ plotList = []
+ for p in peaks:
+ p.append(copy.deepcopy(p[len(p)-1]))
+
+ p[len(p)-1][0] += 1
+
+
+ #print "Adding " + str(p) + " to list"
+ if len(p) >= 0:
+ l = p[len(p)-1]
+ if l[len(l)-1] < 1:
+ with_ = with_.split(" lt")[0] + " lt 9"
+ plotList.append(Gnuplot.Data(p, using="3:1", with_=with_))
+
+
+
+ if len(plotList) > 0 and plot != None:
+ plot(*plotList)
+ time.sleep(0.2)
+
+ #print str(plotList)
+ #for p in peaks:
+ # p = p[0:len(p)-1]
+ return plotList
+
+def main():
+