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
  16 import Gnuplot, Gnuplot.funcutils
  17
  18 gnuplot = Gnuplot.Gnuplot()
  19
  20 def GetData(filename):
  21         input_file = open(filename, "r")
  22         data = []
  23         for line in input_file:
  24                 line = re.sub("#.*", "", line).strip("\r\n ")
  25                 if len(line) == 0:
  26                         continue
  27                 data.append(map(lambda e : float(e), line.split("\t")))
  28         return data
  29
  30 def GetTCS(data):
  31         result = []
  32         n = 0
  33         dI = 0
  34         dE = 0
  35         for i in range(1, len(data)-1):
  36                 dE = data[i+1][1] - data[i][1]
  37                 if (dE != 0):
  38                         n = 0
  39                         dI = 0
  40
  41                 n += 1
  42                 dI += data[i+1][2] - data[i][2]
  43                 if (dE != 0):
  44                         result.append([data[i][1], (dI / (n * dE)) ] ) #/ data[i][2]])
  45         return result
  46
  47 def Plot(*args):
  48         gnuplot.plot(args)
  49
  50 def FitTCS(data):
  51         pass
  52
  53
  54 def main():
  55         if (len(sys.argv) < 2):
  56                 sys.stderr.write(sys.argv[0] + " - Require arguments (filename)\n")
  57                 return 1
  58
  59         tcs = []
  60         gnuplot("set style data lp")
  61         gnuplot("set key outside right")
  62         gnuplot("set title \"Au on Si (95min 3.5A 3-6 e-8mbar)\"")
  63         gnuplot("set xlabel \"E (DAC Counts)\"")
  64         gnuplot("set ylabel \"S(E) (ADC/DAC Counts)\"")
  65         #gnuplot("set term postscript colour")
  66         #gnuplot("set output \"test.eps\"")
  67         for i in range(1, len(sys.argv)):
  68                 tcs.append(GetTCS(GetData(sys.argv[i])))
  69                 if (len(tcs[i-1]) > 0):
  70                         gnuplot.replot(Gnuplot.Data(tcs[i-1], title=sys.argv[i], with_="p"))
  71
  72         # Now average the TCS data
  73         avg = odict.odict()
  74         for t in tcs:
  75                 for p in t:
  76                         if p[0] in avg:
  77                                 avg[p[0]][0] += p[1]
  78                                 avg[p[0]][1] += 1
  79                         else:
  80                                 avg.update({p[0] : [p[1], 1]})
  81
  82         for a in avg.keys():
  83                 avg[a] = float(avg[a][0]) / float(avg[a][1])
  84
  85
  86         gnuplot.replot(Gnuplot.Data(sorted(avg.items(), key = lambda e : e[0]), title="Average", with_="l lw 4"))
  87
  88
  89         print("Press enter to exit")
  90         sys.stdin.readline()
  91
  92         return 0
  93
  94
  95 if __name__ == "__main__":
  96         sys.exit(main())