+ i = 1
+ plotFunc = ShowTCS
+ normalise = False
+ title = ""
+ master_title = ""
+ while i < len(sys.argv):
+ if sys.argv[i] == "--raw":
+ plotFunc = ShowData
+ elif sys.argv[i] == "--tcs":
+ plotFunc = ShowTCS #lambda e : ShowTCS(e, show_peak=False)
+ elif sys.argv[i] == "--output":
+ if i+1 >= len(sys.argv):
+ sys.stderr.write("Need argument for "+sys.argv[i]+" switch\n")
+ sys.exit(1)
+ gnuplot("set term postscript colour")
+ gnuplot("set output \""+sys.argv[i+1]+"\"")
+ i += 1
+ elif sys.argv[i] == "--wxt":
+ gnuplot("set term wxt")
+ elif sys.argv[i] == "--normalise":
+ normalise = True
+ elif sys.argv[i] == "--unnormalise":
+ normalise = False
+ elif sys.argv[i] == "--title":
+ if i+1 >= len(sys.argv):
+ sys.stderr.write("Need argument for "+sys.argv[i]+" switch\n")
+ sys.exit(1)
+ title = sys.argv[i+1]
+ i += 1
+ elif sys.argv[i] == "--master_title":
+ if i+1 >= len(sys.argv):
+ sys.stderr.write("Need argument for "+sys.argv[i]+" switch\n")
+ sys.exit(1)
+ master_title = sys.argv[i+1]
+ i += 1
+ else:
+ plotFunc(sys.argv[i], plot=gnuplot.replot, normalise=normalise, title=title, master_title=master_title)
+
+ i += 1
+
+ print "Done. Press enter to exit, or type name of file to save as."
+ out = sys.stdin.readline().strip("\t\r\n #")
+ if out != "":
+ gnuplot("set term postscript colour")
+ gnuplot("set output \""+out+"\"")
+ gnuplot.replot()
+
+
+def ModelTCS(f, sigma, Emin, Emax, dE):
+ data = []
+ E = Emin
+ while E < Emax:
+ S = (1 - sigma(0))*f(-E) + FuncIntegrate(lambda e : f(e - E) * FuncDerivative(sigma, E, dE), Emin, Emax, dE)
+ data.append([0.00, E, S,0.00])
+ E += dE
+ return data
+
+def IntegrateTCS(data, imin, imax=0, di=1):
+ i = imin
+ if imax == 0:
+ imax = len(data)-1
+ total = 0.0
+
+ while i < imax:
+ total += data[i][2] * (data[i+1][1] - data[i][1])
+ i += di
+ return total
+
+def FuncIntegrate(f, xmin, xmax, dx):
+ x = xmin
+ total = 0.0
+ while x <= xmax:
+ total += f(x) * dx
+ x += dx
+ return total
+
+def FuncDerivative(f, x, dx):
+ return 0.50*(f(x+dx) - f(x-dx))/dx
+
+def FitTCS(data, min_mse=1e-4, max_fail=100, max_adjust=4,divide=10, plot=gnuplot.plot,smooth=0):
+ if type(data) == type(""):
+ d = GetData(data)
+ d = CalibrateData(d)
+ d = MaxNormalise(d)
+ d = Derivative(d)
+ else:
+ d = data
+
+ if smooth != 0:
+ if type(smooth) == type([]):
+ for _ in range(smooth[0]):
+ d = Smooth(d, m=smooth[1])
+ else:
+ d = Smooth(d, m=smooth)
+
+
+
+ plotItems = ShowTCS(d, raw=False,smooth=smooth,plot=None)
+ plotItems.append(None)
+
+ peaks = SmoothPeakFind(d, smooth=5, stop=1, inflection=0)
+ peaks.sort(key = lambda e : e[len(e)-1][1])
+
+ fits = []
+
+ for i in range(0,len(peaks)):
+
+ p = peaks[i]
+ l = p[len(p)-1]
+ if l[len(l)-1] == 0:
+ fits.append([l[3], l[2], 1.0])
+ else:
+ if i-2 >= 0:
+ l = peaks[i-2][len(peaks[i-2])]
+ if l[len(l)-1] == 0:
+ fits.append([l[3], l[2], 1.0])
+ if i+2 <= len(peaks)-1:
+ l = peaks[i+2][len(peaks[i+2])]
+ if l[len(l)-1] == 0:
+ fits.append([l[3], l[2], 1.0])
+
+ for i in range(len(fits)):
+ left = 2.0
+ right = 2.0
+ if i > 0:
+ left = fits[i-1][1] - fits[i][1]
+ if i < len(fits)-1:
+ right = fits[i+1][1] - fits[i][1]
+
+ fits[i][2] = min([abs(0.5*left), abs(0.5*right)])
+
+
+ #print "Fits are " + str(fits)
+ #stdin.readline()
+
+
+
+ def tcs(E):
+ total = 0.0
+ for f in fits:
+ dt = f[0] * gaussian(E - f[1], f[2])
+ #print " Increase total by " + str(dt)
+ total += dt
+ #print "tcs returns " + str(total)
+ return total
+
+ mse = 1
+ old_mse = 1
+ cycle = 0
+ failcount = 0
+
+
+ adjust = 1.0
+
+ iterations = 0
+
+ while failcount < max_fail and mse > min_mse:
+ i = random.randint(0, len(fits)-1)
+ j = random.randint(0, len(fits[i])-1)
+ #while j == 1:
+ # j = random.randint(0, len(fits[i])-1)
+
+ #print "Adjust " + str(i) + ","+str(j) + ": Iteration " + str(iterations) + " mse: " + str(mse)
+
+ old = fits[i][j]
+ old_mse = mse
+
+ fits[i][j] += adjust * (random.random() - 0.50)
+ if i == 2:
+ while fits[i][j] <= 0.0005:
+ fits[i][j] = adjust * (random.random() - 0.50)
+
+
+ model = table(lambda e : [0.00, e, tcs(e), 0.00], 0, 16.8, divide*d[len(d)-1][1]/(len(d)))
+ mse = MeanSquareError(model, d[0::divide])
+ if mse >= old_mse:
+ fits[i][j] = old
+ failcount += 1
+ if failcount > max_fail / 2:
+ if adjust > 1.0/(2.0**max_adjust):
+ adjust /= 2
+ mse = old_mse
+ else:
+ #adjust /= 2.0
+ failcount = 0
+
+
+ iterations += 1