Changes by Chris and project fixes

.gitignore

@@ -0,0 +1,6 @@
+*.class
+*.pyc
+*/__pycache__/*
+nb/PSSS/nbproject/private
+nb/PSSS/build
+nb/PSSS/dist

nb/PSSS/nbproject/private/private.properties

@@ -1,8 +0,0 @@
-compile.on.save=true
-do.depend=false
-do.jar=true
-do.jlink=false
-javac.debug=true
-javadoc.preview=true
-jlink.strip=false
-user.properties.file=/afs/psi.ch/user/g/gobbo_a/.netbeans/12.2/build.properties

nb/PSSS/nbproject/private/private.xml

@@ -1,4 +0,0 @@
-<?xml version="1.0" encoding="UTF-8"?>
-<project-private xmlns="http://www.netbeans.org/ns/project-private/1">
-    <editor-bookmarks xmlns="http://www.netbeans.org/ns/editor-bookmarks/2" lastBookmarkId="0"/>
-</project-private>

script/cpython/psss.py

@@ -0,0 +1,45 @@
+import numpy as np
+from scipy.optimize import curve_fit
+import sys
+
+
+def gaus(x, a, x0, sigma, offset):
+    return offset + a * np.exp(-(x - x0) ** 2 / (2 * sigma ** 2))
+
+#Return [amp, mean_val, sigma, offset]
+def fit_energy(e_from, e_to, steps, num_shots, data):
+    energy_range = np.linspace(e_from, e_to, steps)
+    energy_range_fit = np.linspace(energy_range[0], energy_range[-1], len(energy_range)*10)
+    centre_line_out = data[:,:,int(data.shape[2]/2)].mean(axis=1)
+    try:
+        popt,pcov = curve_fit(gaus,energy_range,centre_line_out,p0=[1,energy_range[np.argmax(centre_line_out)],energy_range.mean()*1e-3,1e3*num_shots])
+    except:
+        raise Exception('Fit failed: spectrum might not be near scan range center \n' + str(sys.exc_info()[1]))
+        #print('Fit failed: spectrum might not be near scan range center')
+        #return None
+    max_ind = np.argmax(centre_line_out)
+    max_photon_energy=energy_range[max_ind]
+    print(max_photon_energy)
+    return popt, centre_line_out
+
+
+#Return [amp, mean_val, sigma, offset]
+def fit_crystal_height(xstal_from, xstal_to, steps, data):
+    xstal_range = np.linspace(xstal_from, xstal_to, steps)
+    projection = data.mean(axis=1).mean(axis=1)
+    offset = np.mean(projection[0:100])
+    signal_centre = xstal_range[np.argmax(projection)]
+    xstal_range_fit = np.linspace(xstal_range[0], xstal_range[-1], len(xstal_range)*10)
+    try:
+        popt,pcov = curve_fit(gaus,xstal_range,projection,p0=[100,signal_centre,-0.2,offset])
+    except:
+        raise Exception('Fit failed: spectrum might not be near scan range center \n' + str(sys.exc_info()[1]))
+        #print('Fit failed: spectrum might not be near scan range center')
+        #return None
+    return popt, projection
+
+
+def get_signal_centre(data, data_range):
+    projection = data.mean(axis=1).mean(axis=1)
+    signal_centre = data_range[np.argmax(projection)]
+    return signal_centre, projection

script/cpython/psss_plot.py

@@ -0,0 +1,32 @@
+import numpy as np
+import matplotlib.pyplot as plt
+
+
+def gaus(x, a, x0, sigma, offset):
+    return offset + a * np.exp(-(x - x0) ** 2 / (2 * sigma ** 2))
+    
+def plot_energy(E_from, E_to, steps, Scan_spec, popt, measured_offset):
+    Energy_range = np.linspace(E_from, E_to, steps)
+    centre_line_out = Scan_spec[:,:,int(Scan_spec.shape[2]/2)].mean(axis=1)
+    Energy_range_fit = np.linspace(Energy_range[0], Energy_range[-1], len(Energy_range)*10)
+
+    
+    plt.figure(figsize=[10,5])
+    plt.subplot(121)
+    plt.title('PSSS scan of set photon energy')
+    plt.pcolormesh(np.arange(0,Scan_spec.shape[2]), Energy_range, Scan_spec.mean(axis=1),cmap='CMRmap')
+    plt.vlines(int(Scan_spec.shape[2]/2), Energy_range[0], Energy_range[-1],linestyles='--', colors='orange')
+    plt.xlim([0,Scan_spec.shape[2]])
+    plt.xlabel('Camera pixel')
+    plt.ylabel('Set PSSS energy [eV] \n SARFE10-PSSS059:ENERGY')
+    
+    plt.subplot(122)
+    plt.title('At camera centre pixel %1i \nCalibrated energy = %.1f [eV]\n Offset from machine = %.1f [eV]'%(int(Scan_spec.shape[2]/2),popt[1],measured_offset))
+    plt.plot(centre_line_out,Energy_range,linewidth = 2, color = 'orange',label ='measured')
+    try:
+        plt.plot(gaus(Energy_range_fit,*popt),Energy_range_fit,'r:',label='fit')
+    except:
+        pass
+    plt.xticks([])
+    plt.legend()
+    plt.grid(True)  

script/cpython/wrapper.py

@@ -0,0 +1,36 @@
+from jeputils import *
+
+RELOAD_CPYTHON = not App.isDetached()
+
+def fit_energy(e_from, e_to, steps, num_shots, data):
+    data = to_array(data, 'd')    
+    dims = [len(data), len(data[0]), len(data[0][0])]
+    data = Convert.flatten(data)
+    arr = to_npa(data, dims)
+    popt, centre_line_out =  call_jep("cpython/psss", "fit_energy", [e_from, e_to, steps, num_shots, arr], reload=RELOAD_CPYTHON)
+    return popt.getData(), centre_line_out.getData()
+
+def fit_crystal_height(xstal_from, xstal_to, steps, data):
+    data = to_array(data, 'd')    
+    dims = [len(data), len(data[0]), len(data[0][0])]
+    data = Convert.flatten(data)
+    arr = to_npa(data, dims)
+    popt, projection =  call_jep("cpython/psss", "fit_crystal_height", [xstal_from, xstal_to, steps, arr], reload=RELOAD_CPYTHON)
+    return popt.getData(), projection.getData()
+
+def get_signal_centre(data, data_range):
+    data = to_array(data, 'd')    
+    dims = [len(data), len(data[0]), len(data[0][0])]
+    data = Convert.flatten(data)   
+    arr = to_npa(data, dims)
+    data_range = to_npa(to_array(data_range, 'd'))
+    signal_centre, projection =  call_jep("cpython/psss", "get_signal_centre", [arr, data_range], reload=RELOAD_CPYTHON)
+    return signal_centre, projection.getData()
+
+def plot_energy(e_from, e_to, steps, data, popt, measured_offset):
+    data = to_array(data, 'd')    
+    dims = [len(data), len(data[0]), len(data[0][0])]
+    data = Convert.flatten(data)
+    arr = to_npa(data, dims)
+    ret =  call_jep("cpython/psss_plot", "plot_energy", [e_from, e_to, steps, arr, popt, measured_offset], reload=RELOAD_CPYTHON)
+    return ret

script/local.py

@@ -2,6 +2,183 @@
 #  Deployment specific global definitions - executed after startup.py
 ###################################################################################################
 
+from mathutils import estimate_peak_indexes, fit_gaussians, create_fit_point_list
+from mathutils import fit_polynomial,fit_gaussian, fit_harmonic, calculate_peaks, fit_gaussian_offset
+from mathutils import PolynomialFunction, Gaussian, HarmonicOscillator, GaussianOffset
+from plotutils import plot_function, plot_data
+
+import java.awt.Color as Color
 run("psss/psss")
 
+###################################################################################################
+# DRY RUN
+###################################################################################################
 
+def set_dry_run(value):
+    global dry_run    
+    dry_run = value
+
+def is_dry_run():
+    if "dry_run" in globals():
+        return True if dry_run else False
+    return False
+          
+
+###################################################################################################
+# Machine utilities
+###################################################################################################
+
+def is_laser_on():
+    return (caget ("SIN-TIMAST-TMA:Beam-Las-Delay-Sel",'d') == 0 )
+
+def is_timing_ok():
+    return caget("SIN-TIMAST-TMA:SOS-COUNT-CHECK") == 0
+
+def get_repetition_rate():
+    return caget("SIN-TIMAST-TMA:Evt-15-Freq-I")
+
+
+###################################################################################################
+# Shortcut to maths utilities
+###################################################################################################
+
+def gfit(ydata, xdata = None):
+    """
+    Gaussian fit
+    """
+    if xdata is None:
+        xdata = frange(0, len(ydata), 1)
+    #ydata = to_list(ydata)
+    #xdata = to_list(xdata)
+    max_y= max(ydata)
+    index_max = ydata.index(max_y)
+    max_x= xdata[index_max]
+    print "Max index:" + str(index_max),
+    print " x:" + str(max_x),
+    print " y:" + str(max_y)
+    gaussians = fit_gaussians(ydata, xdata, [index_max,])
+    (norm, mean, sigma) = gaussians[0]
+    p = plot([ydata],["data"],[xdata], title="Fit" )[0]
+    fitted_gaussian_function = Gaussian(norm, mean, sigma)
+    scale_x = [float(min(xdata)), float(max(xdata)) ]
+    points = max((len(xdata)+1), 100)
+    resolution = (scale_x[1]-scale_x[0]) / points
+    fit_y = []
+    fit_x = frange(scale_x[0],scale_x[1],resolution, True)
+    for x in fit_x:
+        fit_y.append(fitted_gaussian_function.value(x))
+    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)
+        return (None, None, None)
+
+
+def hfit(ydata, xdata = None):
+    """
+    Harmonic fit
+    """
+    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]
+
+    start,end = min(xdata), max(xdata)
+    (amplitude, angular_frequency, phase) = fit_harmonic(ydata, xdata)
+    fitted_harmonic_function = HarmonicOscillator(amplitude, angular_frequency, phase)
+
+    print "amplitude = ", amplitude
+    print "angular frequency = ", angular_frequency
+    print "phase = ", phase
+
+    f = angular_frequency/ (2* math.pi)
+    print "frequency = ", f
+
+    resolution = 4.00 # 1.00
+    fit_y = []
+    for x in frange(start,end,resolution, True):
+        fit_y.append(fitted_harmonic_function.value(x))
+    fit_x = frange(start, end+resolution, resolution)
+
+    p = plot(ydata,"data", xdata, title="HFit")[0]
+    p.addSeries(LinePlotSeries("fit"))
+    p.getSeries(1).setData(fit_x, fit_y)
+
+    #m = (phase + math.pi)/ angular_frequency
+    m = -phase / angular_frequency
+    if (m<start):
+        m+=(1.0/f)
+
+    if start <=m <=end:
+        print "fit = ", m
+        p.addMarker(m, None, "Fit="+str(round(m ,2)), Color.MAGENTA.darker())
+        return (amplitude, angular_frequency, phase, True, m, fit_x, fit_y)
+    else:
+        print "max = ",max_x
+        p.addMarker(max_x, None, "Max="+str(round(max_x ,2)), Color.MAGENTA.darker())
+        return (amplitude, angular_frequency, phase, False, max_x, fit_x, fit_y)
+
+
+
+def plot_gauss_fit(xdata, ydata, gauss_pars=None, p=None, title = "Data"):
+    if gauss_pars is None:
+        gauss_pars= fit_gaussian_offset(ydata, xdata, None)
+    (offset, amp, mean_value, sigma) = gauss_pars
+    print "Gauss plot: ", (offset, amp, mean_value, sigma)
+    fitted_gaussian_function = GaussianOffset(offset, amp, mean_value, abs(sigma))
+    
+    if p is None:
+        p = plot(None, title=title)[0]
+    p.clear()
+    plot_data(p, ydata, title, xdata=xdata,  show_points = True, color=Color.BLUE)
+    fit_range = frange(xdata[0],xdata[-1],float(xdata[1]-xdata[0])/100, True)
+    plot_function(p, fitted_gaussian_function, "Gauss", fit_range, show_points=False, color=Color.RED)
+    p.setLegendVisible(True)
+    p.addMarker(mean_value, None, "Mean=" + str(round(mean_value,2)), Color.LIGHT_GRAY)
+    return p,(amp, mean_value, sigma)
+
+###################################################################################################
+# Tools
+###################################################################################################
+
+def elog(title, message, attachments = [], author = None, category = "Info", domain = "", logbook = "Bernina", encoding=1):
+    """
+    Add entry to ELOG.
+    """
+    if author is None:
+        author = "pshell" #get_context().user.name
+    typ = "pshell"
+    entry = ""
+
+    cmd =           'G_CS_ELOG_add -l "' + logbook       + '" '
+    cmd =     cmd + '-a "Author='        + author        + '" '
+    cmd =     cmd + '-a "Type='          + typ           + '" '
+    cmd =     cmd + '-a "Entry='         + entry         + '" '
+    cmd =     cmd + '-a "Title='         + title         + '" '
+    cmd =     cmd + '-a "Category='      + category      + '" '
+    cmd =     cmd + '-a "Domain='        + domain        + '" '
+    for attachment in attachments:
+        cmd = cmd + '-f "'               + attachment    + '" '
+    cmd =     cmd + '-n '                + str(encoding)
+    cmd =     cmd + ' "'                 + message       + '"'
+    #print cmd
+    #os.system (cmd)
+    #print os.popen(cmd).read()
+    import subprocess
+    proc = subprocess.Popen(cmd, stdout=subprocess.PIPE, shell=True)
+    (out, err) = proc.communicate()
+    if (err is not None) and err!="":
+        raise Exception(err)
+    print out
+    try:
+        return int(out[out.find("ID=") +3 : ])
+    except:
+        print out

script/psss/CrystalHeightScan.py

@@ -10,6 +10,11 @@ if get_exec_pars().source == CommandSource.ui:
     STEPS = 10 #20
     NUM_SHOTS= 10 # 100
     PLOT=None
+# get current camera ROIs and then set to max for scan
+roi_min = psss_roi_min.read()
+roi_max = psss_roi_max.read()
+psss_roi_min.write(1)
+psss_roi_max.write(2000)
 
 p = plot(None, title="Data")[0] if (PLOT is None) else PLOT
 p.clear()
@@ -21,10 +26,10 @@ run("cpython/wrapper")
 
 
 #Setup and functions setup¶
-if not is_dry_run():
-    xstal_height=Channel("SARFE10-PSSS059:MOTOR_Y3.VAL", name="xstal_height")
-else:
-    xstal_height=DummyRegister("xstal_height")
+#if not is_dry_run(): # C.arrell commented out 20.01.21
+xstal_height=Channel("SARFE10-PSSS059:MOTOR_Y3.VAL", name="xstal_height")
+#else:
+#    xstal_height=DummyRegister("xstal_height")
 
 av = create_averager(psss_spectrum_y, NUM_SHOTS, interval=-1, name="spectrum_average")
 av_samples = av.samples
@@ -54,25 +59,12 @@ if not (RANGE_FROM < mean_val < RANGE_TO or RANGE_TO < mean_val < RANGE_FROM):
 xstal_height.write(mean_val)
 xstal_height.close()
 
+# return ROI to inital value
+psss_roi_min.write(roi_min)
+psss_roi_max.write(roi_max)
+
 #Plots
-"""    
-plt.figure(figsize=[10,5])
-plt.subplot(121)
-plt.title('PSSS scan of crystal height')
-plt.pcolormesh(energy_axis, xstal_range, Scan_spec.mean(axis=1),cmap='CMRmap')
-plt.xlim([energy_axis[0],energy_axis[-1]])
-plt.ylim([xstal_range[0], xstal_range[-1]])
-plt.xlabel('PSSS energy axis')
-plt.ylabel('Set crystal position [mm] \n'+PSSS_xstal_height_name[0:-4])
-plt.subplot(122)
-plt.plot(projection,xstal_range,linewidth = 2, color = 'orange',label ='projected signal')
-plt.plot(gaus(xstal_range_fit,*popt),xstal_range_fit,'r:',label='fit')
-plt.ylim([xstal_range[0], xstal_range[-1]])
-plt.title('Signal max at %.3f [mm] (from fit)'%popt[1])
-plt.xticks([])
-plt.legend()
-plt.grid(True)
-"""
+
 p.clear()
 p.setTitle("")
 plot_gauss_fit(xstal_range, projection, gauss_pars=(offset, amp, mean_val, sigma), p=p, title = "Data")

script/psss/EnergyScan.py

@@ -26,15 +26,19 @@ if RANGE_OFF is not None:
 
 run("cpython/wrapper")
 
-
+# get current camera ROIs and then set to max for scan
+roi_min = psss_roi_min.read()
+roi_max = psss_roi_max.read()
+psss_roi_min.write(1)
+psss_roi_max.write(2000)
 
 
 #Scan and take data
 class PSSS_energy(Writable):
     def write(self, value):
-        if not is_dry_run():
-            psss_energy.write(value)
-            exec_cpython("/ioc/modules/qt/PSSS_motion.py", args = ["-m1", "SARFE10-PSSS059"])
+        #if not is_dry_run():
+        psss_energy.write(value)
+        exec_cpython("/ioc/modules/qt/PSSS_motion.py", args = ["-m1", "SARFE10-PSSS059"])
             # python / ioc / modules / qt / PSSS_motion.py - m1 SARFE10 - PSSS059
         time.sleep(1)
         print(value)
@@ -53,10 +57,14 @@ def after_read(record, scan):
       
 r = lscan(en, (av, av_samples), RANGE_FROM, RANGE_TO, STEPS, latency=0.0, after_read = after_read, save=False )
 average, samples, energy_range =  r.getReadable(0), r.getReadable(1), r.getPositions(0)    
+# return ROI to inital value
+psss_roi_min.write(roi_min)
+psss_roi_max.write(roi_max)
 
-[amp, mean_val, sigma, offset],centre_line_out  = fit_energy(RANGE_FROM, RANGE_TO, STEPS+1, NUM_SHOTS, samples)
+[amp, mean_val, sigma, offset],centre_line_out = fit_energy(RANGE_FROM, RANGE_TO, STEPS+1, NUM_SHOTS, samples)
 
 if not (RANGE_FROM < mean_val < RANGE_TO or RANGE_TO < mean_val < RANGE_FROM):
+
     raise Exception ("Invalid fit mean: " + str(mean_val))  
 
     
@@ -67,26 +75,6 @@ print "measured offset", measured_offset
 en.write(mean_val)
 
 
-#Plot
-"""
-plt.figure(figsize=[10,5])
-plt.subplot(121)
-plt.title('PSSS scan of set photon energy')
-plt.pcolormesh(np.arange(0,Scan_spec.shape[2]), Energy_range, Scan_spec.mean(axis=1),cmap='CMRmap')
-plt.vlines(int(Scan_spec.shape[2]/2), Energy_range[0], Energy_range[-1],linestyles='--', colors='orange')
-plt.xlim([0,Scan_spec.shape[2]])
-plt.xlabel('Camera pixel')
-plt.ylabel('Set PSSS energy [eV] \n SARFE10-PSSS059:ENERGY')
-
-plt.subplot(122)
-plt.title('At camera centre pixel %1i \nCalibrated energy = %.1f [eV]\n Offset from machine = %.1f [eV]'%(int(Scan_spec.shape[2]/2),popt[1],measured_offset))
-plt.plot(centre_line_out,Energy_range,linewidth = 2, color = 'orange',label ='measured')
-plt.plot(gaus(Energy_range_fit,*popt),Energy_range_fit,'r:',label='fit')
-plt.xticks([])
-plt.legend()
-plt.grid(True)
-"""
-
 p.clear()
 p.setTitle("")
 plot_gauss_fit(energy_range, centre_line_out, gauss_pars=(offset, amp, mean_val, sigma), p=PLOT, title = "Data")