From: Sam Moore Date: Wed, 10 Oct 2012 23:49:25 +0000 (+0800) Subject: TCS - Improve process.py script X-Git-Url: https://git.ucc.asn.au/?a=commitdiff_plain;h=95832ff21f52524b602dcb863a064873555a1ee9;p=matches%2Fhonours.git TCS - Improve process.py script Going to use this for everything to do with data analysis Abandoning the terrible plot.sh scripts due to fits of rage induced by bash special character escaping --- diff --git a/research/TCS/process.py b/research/TCS/process.py index ff186fb7..31c98d2d 100755 --- a/research/TCS/process.py +++ b/research/TCS/process.py @@ -31,41 +31,85 @@ def GetData(filename): def DoNothing(data): return data -def GetDataSets(directory=".", function=DoNothing): +def AverageAllDataSets(directory=".", function=DoNothing): dirs = {} for f in os.listdir(directory): if os.path.isdir(directory+"/"+str(f)) == True: data_set = [] for datafile in os.listdir(directory+"/"+str(f)): if datafile.split(".")[1] == "dat": - data_set.append(function(map(lambda e : [e[1], e[2]], GetData("./"+str(f)+"/"+str(datafile))))) + data_set.append(GetData(f)) avg = Average(data_set) dirs.update({f : avg}) return dirs +def GetDataSets(directory="."): + data_sets = [] + for f in os.listdir(directory): + if os.path.isdir(directory+"/"+str(f)) == False: + if (len(f.split(".")) > 1 and f.split(".")[1] == "dat"): + d = GetData(directory+"/"+str(f)) + if len(d) > 0: + data_sets.append(d) + return data_sets + -def Derivative(data, a=0, b=1): - result = [] +def Derivative(data, a=1, b=2, sigma=None,step=1): + result = [[]] n = 0 - dI = 0 - dE = 0 - for i in range(1, len(data)-1): - dE = data[i+1][a] - data[i][a] - if (dE != 0): - n = 0 - dI = 0 - - n += 1 - dI += data[i+1][b] - data[i][b] - if (dE != 0): - result.append([data[i][a], (dI / (n * dE)) ] ) #/ data[i][2]]) - return result + dI = [0,0] + dE = [0,0] + + for i in range(0, len(data),step): + result[len(result)-1] = [d for d in data[i]] + if (i >= step): + dE[0] = data[i][a] - data[i-step][a] + dI[0] = data[i][b] - data[i-step][b] + else: + dI[0] = None -def MaxNormalise(data, u=1): + if (i < len(data)-step): + dE[1] = data[i+step][a] - data[i][a] + dI[1] = data[i+step][b] - data[i][b] + else: + dI[1] = None - + if sigma != None: + #print str(data[i]) + " ["+str(sigma)+"] = " + str(data[i][int(abs(sigma))]) + if sigma < 0: + if dI[0] != None: dI[0] -= 0.5*data[i][int(abs(sigma))] + if dI[1] != None: dI[1] -= 0.5*data[i][int(abs(sigma))] + else: + if dI[0] != None: dI[0] += 0.5*data[i][int(abs(sigma))] + if dI[1] != None: dI[1] += 0.5*data[i][int(abs(sigma))] + + deltaE = 0.0 + deltaI = 0.0 + count = 0 + if dI[0] != None: + deltaE += dE[0] + deltaI += dI[0] + count += 1 + if dI[1] != None: + deltaE += dE[1] + deltaI += dI[1] + count += 1 + + deltaI /= float(count) + deltaE /= float(count) + + + if (deltaE != 0): + result[len(result)-1][b] = (deltaI / deltaE) + else: + result[len(result)-1][b] = 0.0 + result.append([]) + + return result[0:len(result)-1] + +def MaxNormalise(data, u=1): result = copy.deepcopy(data) if (len(data) <= 0): return result @@ -80,21 +124,19 @@ def Average(data_sets, u=1): avg = odict.odict() for t in data_sets: for p in t: - if p[0] in avg: - avg[p[0]][0] += p[u] - avg[p[0]][1] += 1 + if p[u] in avg: + #print "Already have " + str(p[u]) + avg[p[u]][1] += 1 + for i in range(0, len(p)): + avg[p[u]][0][i] += p[i] else: - avg.update({p[0] : [p[u], 1]}) + #print "Create key for " + str(p[u]) + avg.update({p[u] : [p, 1]}) for a in avg.keys(): - avg[a] = float(avg[a][0]) / float(avg[a][1]) - return sorted(avg.items(), key = lambda e : e[0]) - -def Plot(*args): - gnuplot.plot(args) - -def FitTCS(data): - pass + for i in range(0, len(avg[a][0])): + avg[a][0][i] /= float(avg[a][1]) + return map(lambda e : e[1][0], sorted(avg.items(), key = lambda e : e[0])) def FullWidthAtHalfMax(data, u=1): maxval = max(data, key = lambda e : e[u]) @@ -119,8 +161,92 @@ def FullWidthAtHalfMax(data, u=1): def SaveData(filename, data): out = open(filename, "w", 0) for a in data: - out.write(str(a[0]) + "\t" + str(a[1]) + "\n") + for i in range(0, len(a)): + out.write(str(a[i])) + if (i < len(a) - 1): + out.write("\t") + out.write("\n") + +def CalibrateData(data, ammeter_scale=1e-6): + 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, calibrate=True, normalise=False, show_error=False, plot=gnuplot.plot,w="lp", step=1): + if type(filename) == type(""): + data = GetData(filename) + else: + data = filename + filename = "data" + + 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)\"") + + gnuplot("set title \"S(E)\"") + gnuplot("set xlabel \"U ("+str(units[0])+")\"") + + + d = Derivative(data, 1, 2, step=step) + + plot(Gnuplot.Data(d, using="2:3", with_=w,title="S(E) : " + str(filename))) + if (show_error): + error1 = Derivative(data, 1, 2, -3,step=step) + error2 = Derivative(data, 1, 2, +3,step=step) + gnuplot.replot(Gnuplot.Data(error1, using="2:3", with_=w,title="Error : Low bound")) + gnuplot.replot(Gnuplot.Data(error2, using="2:3", with_=w, title="Error : Upper bound")) + + return data + +def ShowData(filename,calibrate=True, normalise=False, show_error=False, plot=gnuplot.plot,w="lp", step=1): + if type(filename) == type(""): + data = GetData(filename) + else: + data = filename + filename = "data" + + if len(data) <= 0: + return + 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)\"") + + gnuplot("set title \"S(E)\"") + gnuplot("set xlabel \"U ("+str(units[0])+")\"") + + + #d = Derivative(data, 1, 2, step=step) + + plot(Gnuplot.Data(data, using="2:3", with_=w,title="S(E) : " + str(filename))) + 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]) + gnuplot.replot(Gnuplot.Data(error1, using="2:3", with_=w,title="Error : Low bound")) + gnuplot.replot(Gnuplot.Data(error2, using="2:3", with_=w, title="Error : Upper bound")) + + return data def main(): if (len(sys.argv) < 2): @@ -136,7 +262,7 @@ def main(): #gnuplot("set term postscript colour") #gnuplot("set output \"test.eps\"") for i in range(1, len(sys.argv)): - tcs.append(DoNothing(map(lambda e : [e[1], e[2]], GetData(sys.argv[i])))) + tcs.append(Derivative(GetData(sys.argv[i]), 1, 2)) #tcs.append(GetTCS(GetData(sys.argv[i]))) if (len(tcs[i-1]) > 0): gnuplot.replot(Gnuplot.Data(tcs[i-1], title=sys.argv[i], with_="lp")) @@ -146,7 +272,9 @@ def main(): avg = Average(tcs) - #gnuplot.replot(Gnuplot.Data(avg, title="Average", with_="l lw 2")) + for a in avg: + sys.stdout.write(str(a[0]) + "\t" + str(a[1]) + "\t" + str(a[1]) + "\n") + gnuplot.replot(Gnuplot.Data(avg, title="Average", with_="l lw 2")) sys.stdout.write("Save averaged data as (blank for no save): ") filename = sys.stdin.readline().strip(" \r\n\t")