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 GetDataSets(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(function(map(lambda e : [e[1], e[2]], GetData("./"+str(f)+"/"+str(datafile)))))
  42
  43                         avg = Average(data_set)
  44                         dirs.update({f : avg})
  45         return dirs
  46
  47
  48
  49 def Derivative(data, a=0, b=1):
  50         result = []
  51         n = 0
  52         dI = 0
  53         dE = 0
  54         for i in range(1, len(data)-1):
  55                 dE = data[i+1][a] - data[i][a]
  56                 if (dE != 0):
  57                         n = 0
  58                         dI = 0
  59
  60                 n += 1
  61                 dI += data[i+1][b] - data[i][b]
  62                 if (dE != 0):
  63                         result.append([data[i][a], (dI / (n * dE)) ] ) #/ data[i][2]])
  64         return result
  65
  66 def MaxNormalise(data, u=1):
  67
  68
  69         result = copy.deepcopy(data)
  70         if (len(data) <= 0):
  71                 return result
  72         maxval = max(data, key = lambda e : e[u])[u]
  73
  74         for d in result:
  75                 d[u] = d[u] / maxval
  76
  77         return result
  78
  79 def Average(data_sets, u=1):
  80         avg = odict.odict()
  81         for t in data_sets:
  82                 for p in t:
  83                         if p[0] in avg:
  84                                 avg[p[0]][0] += p[u]
  85                                 avg[p[0]][1] += 1
  86                         else:
  87                                 avg.update({p[0] : [p[u], 1]})
  88
  89         for a in avg.keys():
  90                 avg[a] = float(avg[a][0]) / float(avg[a][1])
  91         return sorted(avg.items(), key = lambda e : e[0])
  92
  93 def Plot(*args):
  94         gnuplot.plot(args)
  95
  96 def FitTCS(data):
  97         pass
  98
  99 def FullWidthAtHalfMax(data, u=1):
 100         maxval = max(data, key = lambda e : e[u])
 101         peak = data.index(maxval)
 102         maxval = maxval[0]
 103         lhs = None
 104         rhs = None
 105         for i in range(1, len(data)/2):
 106                 if lhs == None:
 107                         if (peak-i > 0 and data[peak-i] < 0.50*maxval):
 108                                 lhs = data[peak-i][u]
 109                 if rhs == None:
 110                         if (peak+i < len(data) and data[peak+i] < 0.50*maxval):
 111                                 rhs = peak+i
 112                 if lhs != None and rhs != None:
 113                         break
 114         if rhs == None or lhs == None:
 115                 return abs(data[len(data)-1][0] - data[0][0])
 116         else:
 117                 return abs(rhs - lhs)
 118
 119 def SaveData(filename, data):
 120         out = open(filename, "w", 0)
 121         for a in data:
 122                 out.write(str(a[0]) + "\t" + str(a[1]) + "\n")
 123
 124
 125 def main():
 126         if (len(sys.argv) < 2):
 127                 sys.stderr.write(sys.argv[0] + " - Require arguments (filename)\n")
 128                 return 1
 129
 130         tcs = []
 131         gnuplot("set style data lp")
 132         gnuplot("set key outside right")
 133         #gnuplot("set title \"Au on Si (50min 3.5A 3-6 e-8mbar)\"")
 134         #gnuplot("set xlabel \"E (DAC Counts)\"")
 135         #gnuplot("set ylabel \"S(E) (ADC/DAC Counts)\"")
 136         #gnuplot("set term postscript colour")
 137         #gnuplot("set output \"test.eps\"")
 138         for i in range(1, len(sys.argv)):
 139                 tcs.append(DoNothing(map(lambda e : [e[1], e[2]], GetData(sys.argv[i]))))
 140                 #tcs.append(GetTCS(GetData(sys.argv[i])))
 141                 if (len(tcs[i-1]) > 0):
 142                         gnuplot.replot(Gnuplot.Data(tcs[i-1], title=sys.argv[i], with_="lp"))
 143
 144         # Now average the data
 145
 146
 147
 148         avg = Average(tcs)
 149         #gnuplot.replot(Gnuplot.Data(avg, title="Average", with_="l lw 2"))
 150
 151         sys.stdout.write("Save averaged data as (blank for no save): ")
 152         filename = sys.stdin.readline().strip(" \r\n\t")
 153         if (filename != ""):
 154                 SaveData(filename, avg)
 155
 156         return 0
 157
 158
 159 if __name__ == "__main__":
 160         sys.exit(main())