Closedown

script/local.py

@@ -16,3 +16,45 @@ def get_energy():
     print 'Energy [keV]:'+ str(e) + '   Wavelength [A]:' + str(12.39842/e)
     return e
 
+
+#Fitting
+from mathutils import estimate_peak_indexes, fit_gaussians, create_fit_point_list, Gaussian
+import java.awt.Color as Color
+
+def fit(ydata, xdata = None):
+    if xdata is None:
+        xdata = frange(0, len(ydata), 1)
+    max_y= max(ydata)    
+    index_max = ydata.index(max_y)    
+    max_x= xdata[index_max]
+    print "Max indexs: " + str(index_max),
+    print " Max x: " + str(max_x),
+    print " Max y: " + str(max_y)
+    
+    gaussians = fit_gaussians(ydata, xdata, [index_max,])
+    (norm, mean, sigma) = gaussians[0]
+
+    p = plot([ydata],["data"],[xdata], context="Fit" )[0] 
+    fitted_gaussian_function = Gaussian(norm, mean, sigma)
+
+    scale_x = [float(min(xdata)), float(max(xdata)) ]
+    resolution = (scale_x[1]-scale_x[0]) / 100
+    print resolution
+    fit_y = []
+    for x in frange(scale_x[0],scale_x[1],resolution, True):
+        fit_y.append(fitted_gaussian_function.value(x))
+    fit_x = frange(scale_x[0], scale_x[1]+resolution, resolution)
+    #plots = plot([ydata, fit_y] ,["data", "fit"], xdata = [xdata,fit_x], context="Fit")             
+    p.addSeries(LinePlotSeries("fit"))
+    p.getSeries(1).setData(fit_x, fit_y)    
+    
+    if abs(mean - xdata[index_max]) < ((scale_x[0] + scale_x[1])/2):
+        print "Mean -> " +  str(mean)
+        p.addMarker(mean, None, "Mean="+str(round(norm,2)), Color.MAGENTA.darker())
+        return (norm, mean, sigma)
+    else:
+        p.addMarker(max_x, None, "Max="+str(round(max_x,2)), Color.Gray)
+        print "Invalid gaussian fit: " +  str(mean)
+    
+
+#fit([1,2,3, 3,3,3,4,5,6,12,33,23,15,3,2,1, 1, 1, 1, 1])