research/TCS/process.py

   1 #!/usr/bin/python -u
   2
   3 #
   4 # @file process.py
   5 # @purpose Process TCS data
   6 #               Takes S(E) = dI/dE
   7 # @author Sam Moore
   8 # @date August 2012
   9 #
  10
  11 import sys
  12 import os
  13 import re # Regular expressions - for removing comments
  14 import odict #ordered dictionary
  15 import copy
  16
  17 import Gnuplot, Gnuplot.funcutils
  18 import string
  19
  20 gnuplot = Gnuplot.Gnuplot()
  21
  22 def Reset():
  23         gnuplot = Gnuplot.Gnuplot()
  24
  25 def FindDataFiles(directory=".", depth=1, result=None):
  26         if result == None:
  27                 result = []
  28
  29         for f in os.listdir(directory):
  30                 if os.path.isdir(directory+"/"+str(f)):
  31                         if depth > 1:
  32                                 result += FindDataFiles(directory+"/"+str(f), depth-1, result)
  33                         continue
  34                 s = f.split(".")
  35                 if (len(s) == 2 and s[1] == "dat"):
  36                         result.append(directory+"/"+str(f))
  37
  38         return result
  39
  40 def BaseName(f):
  41         a = f.split("/")
  42         return a[len(a)-1]
  43
  44
  45 def DirectoryName(f, start=0,back=1):
  46         a = f.split("/")
  47         return string.join(a[start:(len(a)-back)], "/")
  48
  49 def GetData(filename):
  50         input_file = open(filename, "r")
  51         data = []
  52         for line in input_file:
  53                 line = re.sub("#.*", "", line).strip("\r\n ")
  54                 if len(line) == 0:
  55                         continue
  56                 data.append(map(lambda e : float(e), line.split("\t")))
  57         return data
  58
  59 def DoNothing(data):
  60         return data
  61
  62 def AverageAllDataSets(directory=".", function=DoNothing):
  63         dirs = {}
  64         for f in os.listdir(directory):
  65                 if os.path.isdir(directory+"/"+str(f)) == True:
  66                         data_set = []
  67                         for datafile in os.listdir(directory+"/"+str(f)):
  68                                 if datafile.split(".")[1] == "dat":
  69                                         data_set.append(GetData(f))
  70
  71                         avg = Average(data_set)
  72                         dirs.update({f : avg})
  73         return dirs
  74
  75 def GetDataSets(directory="."):
  76         data_sets = []
  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))
  81                                 if len(d) > 0:
  82                                         data_sets.append(d)
  83         return data_sets
  84
  85
  86
  87 def Derivative(data, a=1, b=2, sigma=None,step=1):
  88         result = [[]]
  89         n = 0
  90         dI = [0,0]
  91         dE = [0,0]
  92
  93         for i in range(0, len(data),step):
  94                 result[len(result)-1] = [d for d in data[i]]
  95                 if (i >= step):
  96                         dE[0] = data[i][a] - data[i-step][a]
  97                         dI[0] = data[i][b] - data[i-step][b]
  98                 else:
  99                         dI[0] = None
 100
 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]
 104                 else:
 105                         dI[1] = None
 106
 107                 if sigma != None:
 108                         #print str(data[i]) + " ["+str(sigma)+"] = " + str(data[i][int(abs(sigma))])
 109                         if sigma < 0:
 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))]
 112                         else:
 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))]
 115
 116                 deltaE = 0.0
 117                 deltaI = 0.0
 118                 count = 0
 119                 if dI[0] != None:
 120                         deltaE += dE[0]
 121                         deltaI += dI[0]
 122                         count += 1
 123                 if dI[1] != None:
 124                         deltaE += dE[1]
 125                         deltaI += dI[1]
 126                         count += 1
 127
 128                 deltaI /= float(count)
 129                 deltaE /= float(count)
 130
 131
 132                 if (deltaE != 0):
 133                         result[len(result)-1][b] = (deltaI / deltaE)
 134                 else:
 135                         result[len(result)-1][b] = 0.0
 136                 result.append([])
 137
 138         return result[0:len(result)-1]
 139
 140 def MaxNormalise(data, u=2):
 141         result = copy.deepcopy(data)
 142         if (len(data) <= 0):
 143                 return result
 144         maxval = max(data, key = lambda e : e[u])[u]
 145
 146         for d in result:
 147                 d[u] = d[u] / maxval
 148
 149         return result
 150
 151 def Average(data_sets, u=1):
 152         avg = odict.odict()
 153         for t in data_sets:
 154                 for p in t:
 155                         if p[u] in avg:
 156                                 #print "Already have " + str(p[u])
 157                                 avg[p[u]][1] += 1
 158                                 for i in range(0, len(p)):
 159                                         avg[p[u]][0][i] += p[i]
 160                         else:
 161                                 #print "Create key for " + str(p[u])
 162                                 avg.update({p[u] : [p, 1]})
 163
 164         for a in avg.keys():
 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]))
 168
 169 def FullWidthAtHalfMax(data, u=1):
 170         maxval = max(data, key = lambda e : e[u])
 171         peak = data.index(maxval)
 172         maxval = maxval[0]
 173         lhs = None
 174         rhs = None
 175         for i in range(1, len(data)/2):
 176                 if lhs == None:
 177                         if (peak-i > 0 and data[peak-i] < 0.50*maxval):
 178                                 lhs = data[peak-i][u]
 179                 if rhs == None:
 180                         if (peak+i < len(data) and data[peak+i] < 0.50*maxval):
 181                                 rhs = peak+i
 182                 if lhs != None and rhs != None:
 183                         break
 184         if rhs == None or lhs == None:
 185                 return abs(data[len(data)-1][0] - data[0][0])
 186         else:
 187                 return abs(rhs - lhs)
 188
 189 def SaveData(filename, data):
 190         out = open(filename, "w", 0)
 191         for a in data:
 192                 for i in range(0, len(a)):
 193                         out.write(str(a[i]))
 194                         if (i < len(a) - 1):
 195                                 out.write("\t")
 196                 out.write("\n")
 197
 198 def AverageAllData(directory, save=None)
 199         data_sets = []
 200         if save == None: save = directory+"/average.dat":
 201         for d in FindDataFiles(directory):
 202                 data_sets.append(GetData(d))
 203
 204         a = Average(data_sets)
 205         SaveData(save, a)
 206         return a
 207
 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
 214         return data
 215
 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)
 219         else:
 220                 data = filename
 221                 filename = "tcs data"
 222
 223         if (title == ""):
 224                 title = BaseName(filename)
 225
 226         if (len(data) <= 0):
 227                 return data
 228
 229         if calibrate:
 230                 data = CalibrateData(data)
 231                 units = ["V", "uA / V"]
 232         else:
 233                 units = ["DAC counts", "ADC counts / DAC counts"]
 234
 235         if not normalise:
 236                 gnuplot("set ylabel \"dI(E)/dE ("+str(units[1])+")\"")
 237         else:
 238                 data = MaxNormalise(data)
 239                 gnuplot("set ylabel \"dI(E)/dE (normalised)\"")
 240
 241         if (output != None and type(output) == type("")):
 242                 gnuplot("set term png size 640,480")
 243                 gnuplot("set output \""+str(output)+"\"")
 244
 245         gnuplot("set title \"Total Current Spectrum S(E)\"")
 246         gnuplot("set xlabel \"U ("+str(units[0])+")\"")
 247
 248
 249         d = Derivative(data, 1, 2, step=step)
 250
 251         plot(Gnuplot.Data(d, using="2:3", with_=with_,title=title))
 252         if (show_error):
 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"))
 257
 258         if (output != None and type(output) == type("")):
 259                 gnuplot("set term wxt")
 260         return data
 261
 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)
 265         else:
 266                 data = filename
 267                 filename = "raw data"
 268
 269         if (title == ""):
 270                 title = BaseName(filename)
 271
 272         if len(data) <= 0:
 273                 return data
 274         if calibrate:
 275                 data = CalibrateData(data)
 276                 units = ["V", "uA"]
 277         else:
 278                 units = ["DAC counts", "ADC counts"]
 279
 280         if not normalise:
 281                 gnuplot("set ylabel \"I(E) ("+str(units[1])+")\"")
 282         else:
 283                 data = MaxNormalise(data)
 284                 gnuplot("set ylabel \"I(E) (normalised)\"")
 285
 286         if (output != None and type(output) == type("")):
 287                 gnuplot("set term png size 640,480")
 288                 gnuplot("set output \""+str(output)+"\"")
 289
 290         gnuplot("set title \"Sample Current I(E)\"")
 291         gnuplot("set xlabel \"U ("+str(units[0])+")\"")
 292
 293
 294         #d = Derivative(data, 1, 2, step=step)
 295
 296         plot(Gnuplot.Data(data, using="2:3", with_=with_,title=title))
 297         if (show_error):
 298                 error1 = copy.deepcopy(data)
 299                 error2 = copy.deepcopy(data)
 300                 for i in range(len(data)):
 301                         #print str(data[i])
 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"))
 306
 307         if (output != None and type(output) == type("")):
 308                 gnuplot("set term wxt")
 309         return data
 310
 311 def main():
 312         return 0
 313         if (len(sys.argv) < 2):
 314                 sys.stderr.write(sys.argv[0] + " - Require arguments (filename)\n")
 315                 return 1
 316
 317         tcs = []
 318         gnuplot("set style data lp")
 319         gnuplot("set key outside right")
 320         #gnuplot("set title \"Au on Si (50min 3.5A 3-6 e-8mbar)\"")
 321         #gnuplot("set xlabel \"E (DAC Counts)\"")
 322         #gnuplot("set ylabel \"S(E) (ADC/DAC Counts)\"")
 323         #gnuplot("set term postscript colour")
 324         #gnuplot("set output \"test.eps\"")
 325         for i in range(1, len(sys.argv)):
 326                 if (len(tcs[i-1]) > 0):
 327                         gnuplot.replot(Gnuplot.Data(tcs[i-1], title=sys.argv[i], with_="lp"))
 328
 329         # Now average the data
 330
 331
 332
 333         avg = Average(tcs)
 334         for a in avg:
 335                 sys.stdout.write(str(a[0]) + "\t" + str(a[1]) + "\t" + str(a[1]) + "\n")
 336         gnuplot.replot(Gnuplot.Data(avg, title="Average", with_="l lw 2"))
 337
 338         sys.stdout.write("Save averaged data as (blank for no save): ")
 339         filename = sys.stdin.readline().strip(" \r\n\t")
 340         if (filename != ""):
 341                 SaveData(filename, avg)
 342
 343         return 0
 344
 345
 346 if __name__ == "__main__":
 347         sys.exit(main())