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
  19 gnuplot = Gnuplot.Gnuplot()
  20
  21 def GetData(filename):
  22         input_file = open(filename, "r")
  23         data = []
  24         for line in input_file:
  25                 line = re.sub("#.*", "", line).strip("\r\n ")
  26                 if len(line) == 0:
  27                         continue
  28                 data.append(map(lambda e : float(e), line.split("\t")))
  29         return data
  30
  31 def DoNothing(data):
  32         return data
  33
  34 def AverageAllDataSets(directory=".", function=DoNothing):
  35         dirs = {}
  36         for f in os.listdir(directory):
  37                 if os.path.isdir(directory+"/"+str(f)) == True:
  38                         data_set = []
  39                         for datafile in os.listdir(directory+"/"+str(f)):
  40                                 if datafile.split(".")[1] == "dat":
  41                                         data_set.append(GetData(f))
  42
  43                         avg = Average(data_set)
  44                         dirs.update({f : avg})
  45         return dirs
  46
  47 def GetDataSets(directory="."):
  48         data_sets = []
  49         for f in os.listdir(directory):
  50                 if os.path.isdir(directory+"/"+str(f)) == False:
  51                         if (len(f.split(".")) > 1 and f.split(".")[1] == "dat"):
  52                                 d = GetData(directory+"/"+str(f))
  53                                 if len(d) > 0:
  54                                         data_sets.append(d)
  55         return data_sets
  56
  57
  58
  59 def Derivative(data, a=1, b=2, sigma=None,step=1):
  60         result = [[]]
  61         n = 0
  62         dI = [0,0]
  63         dE = [0,0]
  64
  65         for i in range(0, len(data),step):
  66                 result[len(result)-1] = [d for d in data[i]]
  67                 if (i >= step):
  68                         dE[0] = data[i][a] - data[i-step][a]
  69                         dI[0] = data[i][b] - data[i-step][b]
  70                 else:
  71                         dI[0] = None
  72
  73                 if (i < len(data)-step):
  74                         dE[1] = data[i+step][a] - data[i][a]
  75                         dI[1] = data[i+step][b] - data[i][b]
  76                 else:
  77                         dI[1] = None
  78
  79                 if sigma != None:
  80                         #print str(data[i]) + " ["+str(sigma)+"] = " + str(data[i][int(abs(sigma))])
  81                         if sigma < 0:
  82                                 if dI[0] != None: dI[0] -= 0.5*data[i][int(abs(sigma))]
  83                                 if dI[1] != None: dI[1] -= 0.5*data[i][int(abs(sigma))]
  84                         else:
  85                                 if dI[0] != None: dI[0] += 0.5*data[i][int(abs(sigma))]
  86                                 if dI[1] != None: dI[1] += 0.5*data[i][int(abs(sigma))]
  87
  88                 deltaE = 0.0
  89                 deltaI = 0.0
  90                 count = 0
  91                 if dI[0] != None:
  92                         deltaE += dE[0]
  93                         deltaI += dI[0]
  94                         count += 1
  95                 if dI[1] != None:
  96                         deltaE += dE[1]
  97                         deltaI += dI[1]
  98                         count += 1
  99
 100                 deltaI /= float(count)
 101                 deltaE /= float(count)
 102
 103
 104                 if (deltaE != 0):
 105                         result[len(result)-1][b] = (deltaI / deltaE)
 106                 else:
 107                         result[len(result)-1][b] = 0.0
 108                 result.append([])
 109
 110         return result[0:len(result)-1]
 111
 112 def MaxNormalise(data, u=1):
 113         result = copy.deepcopy(data)
 114         if (len(data) <= 0):
 115                 return result
 116         maxval = max(data, key = lambda e : e[u])[u]
 117
 118         for d in result:
 119                 d[u] = d[u] / maxval
 120
 121         return result
 122
 123 def Average(data_sets, u=1):
 124         avg = odict.odict()
 125         for t in data_sets:
 126                 for p in t:
 127                         if p[u] in avg:
 128                                 #print "Already have " + str(p[u])
 129                                 avg[p[u]][1] += 1
 130                                 for i in range(0, len(p)):
 131                                         avg[p[u]][0][i] += p[i]
 132                         else:
 133                                 #print "Create key for " + str(p[u])
 134                                 avg.update({p[u] : [p, 1]})
 135
 136         for a in avg.keys():
 137                 for i in range(0, len(avg[a][0])):
 138                         avg[a][0][i] /= float(avg[a][1])
 139         return map(lambda e : e[1][0], sorted(avg.items(), key = lambda e : e[0]))
 140
 141 def FullWidthAtHalfMax(data, u=1):
 142         maxval = max(data, key = lambda e : e[u])
 143         peak = data.index(maxval)
 144         maxval = maxval[0]
 145         lhs = None
 146         rhs = None
 147         for i in range(1, len(data)/2):
 148                 if lhs == None:
 149                         if (peak-i > 0 and data[peak-i] < 0.50*maxval):
 150                                 lhs = data[peak-i][u]
 151                 if rhs == None:
 152                         if (peak+i < len(data) and data[peak+i] < 0.50*maxval):
 153                                 rhs = peak+i
 154                 if lhs != None and rhs != None:
 155                         break
 156         if rhs == None or lhs == None:
 157                 return abs(data[len(data)-1][0] - data[0][0])
 158         else:
 159                 return abs(rhs - lhs)
 160
 161 def SaveData(filename, data):
 162         out = open(filename, "w", 0)
 163         for a in data:
 164                 for i in range(0, len(a)):
 165                         out.write(str(a[i]))
 166                         if (i < len(a) - 1):
 167                                 out.write("\t")
 168                         out.write("\n")
 169
 170 def CalibrateData(data, ammeter_scale=1e-6):
 171         for i in range(0, len(data)):
 172                 data[i][1] = 16.8 * float(data[i][1]) / 4000.0
 173                 data[i][2] = ammeter_scale * 0.170 * float(data[i][2]) / 268.0
 174                 data[i][3] = ammeter_scale * 0.170 * float(data[i][3]) / 268.0
 175         return data
 176
 177 def ShowTCS(filename, calibrate=True, normalise=False, show_error=False, plot=gnuplot.plot,w="lp", step=1):
 178         if type(filename) == type(""):
 179                 data = GetData(filename)
 180         else:
 181                 data = filename
 182                 filename = "data"
 183
 184         if calibrate:
 185                 data = CalibrateData(data)
 186                 units = ["V", "uA / V"]
 187         else:
 188                 units = ["DAC counts", "ADC counts / DAC counts"]
 189
 190         if not normalise:
 191                 gnuplot("set ylabel \"dI(E)/dE ("+str(units[1])+")\"")
 192         else:
 193                 data = MaxNormalise(data)
 194                 gnuplot("set ylabel \"dI(E)/dE (normalised)\"")
 195
 196         gnuplot("set title \"S(E)\"")
 197         gnuplot("set xlabel \"U ("+str(units[0])+")\"")
 198
 199
 200         d = Derivative(data, 1, 2, step=step)
 201
 202         plot(Gnuplot.Data(d, using="2:3", with_=w,title="S(E) : " + str(filename)))
 203         if (show_error):
 204                 error1 = Derivative(data, 1, 2, -3,step=step)
 205                 error2 = Derivative(data, 1, 2, +3,step=step)
 206                 gnuplot.replot(Gnuplot.Data(error1, using="2:3", with_=w,title="Error : Low bound"))
 207                 gnuplot.replot(Gnuplot.Data(error2, using="2:3", with_=w, title="Error : Upper bound"))
 208
 209         return data
 210
 211 def ShowData(filename,calibrate=True, normalise=False, show_error=False, plot=gnuplot.plot,w="lp", step=1):
 212         if type(filename) == type(""):
 213                 data = GetData(filename)
 214         else:
 215                 data = filename
 216                 filename = "data"
 217
 218         if len(data) <= 0:
 219                 return
 220         if calibrate:
 221                 data = CalibrateData(data)
 222                 units = ["V", "uA"]
 223         else:
 224                 units = ["DAC counts", "ADC counts"]
 225
 226         if not normalise:
 227                 gnuplot("set ylabel \"I(E) ("+str(units[1])+")\"")
 228         else:
 229                 data = MaxNormalise(data)
 230                 gnuplot("set ylabel \"I(E) (normalised)\"")
 231
 232         gnuplot("set title \"S(E)\"")
 233         gnuplot("set xlabel \"U ("+str(units[0])+")\"")
 234
 235
 236         #d = Derivative(data, 1, 2, step=step)
 237
 238         plot(Gnuplot.Data(data, using="2:3", with_=w,title="S(E) : " + str(filename)))
 239         if (show_error):
 240                 error1 = copy.deepcopy(data)
 241                 error2 = copy.deepcopy(data)
 242                 for i in range(len(data)):
 243                         #print str(data[i])
 244                         error1[i][2] -= 0.50*float(data[i][3])
 245                         error2[i][2] += 0.50*float(data[i][3])
 246                 gnuplot.replot(Gnuplot.Data(error1, using="2:3", with_=w,title="Error : Low bound"))
 247                 gnuplot.replot(Gnuplot.Data(error2, using="2:3", with_=w, title="Error : Upper bound"))
 248
 249         return data
 250
 251 def main():
 252         if (len(sys.argv) < 2):
 253                 sys.stderr.write(sys.argv[0] + " - Require arguments (filename)\n")
 254                 return 1
 255
 256         tcs = []
 257         gnuplot("set style data lp")
 258         gnuplot("set key outside right")
 259         #gnuplot("set title \"Au on Si (50min 3.5A 3-6 e-8mbar)\"")
 260         #gnuplot("set xlabel \"E (DAC Counts)\"")
 261         #gnuplot("set ylabel \"S(E) (ADC/DAC Counts)\"")
 262         #gnuplot("set term postscript colour")
 263         #gnuplot("set output \"test.eps\"")
 264         for i in range(1, len(sys.argv)):
 265 <<<<<<< HEAD
 266                 tcs.append(Derivative(map(lambda e : [e[1], e[2]], GetData(sys.argv[i]))))
 267 =======
 268                 tcs.append(Derivative(GetData(sys.argv[i]), 1, 2))
 269 >>>>>>> 95832ff21f52524b602dcb863a064873555a1ee9
 270                 #tcs.append(GetTCS(GetData(sys.argv[i])))
 271                 if (len(tcs[i-1]) > 0):
 272                         gnuplot.replot(Gnuplot.Data(tcs[i-1], title=sys.argv[i], with_="lp"))
 273
 274         # Now average the data
 275
 276
 277
 278         avg = Average(tcs)
 279         for a in avg:
 280                 sys.stdout.write(str(a[0]) + "\t" + str(a[1]) + "\t" + str(a[1]) + "\n")
 281         gnuplot.replot(Gnuplot.Data(avg, title="Average", with_="l lw 2"))
 282
 283         sys.stdout.write("Save averaged data as (blank for no save): ")
 284         filename = sys.stdin.readline().strip(" \r\n\t")
 285         if (filename != ""):
 286                 SaveData(filename, avg)
 287
 288         return 0
 289
 290
 291 if __name__ == "__main__":
 292         sys.exit(main())