From 33e2f84f5ec5d33fbb530d9fe8e5fe03c7321684 Mon Sep 17 00:00:00 2001 From: gobbo_a Date: Thu, 4 May 2017 12:42:06 +0200 Subject: [PATCH] Startup --- devices/CurrentCamera.properties | 20 +- script/Diagnostics/WireScan.py | 1 + script/Diagnostics/sig_process.py | 223 ++++++++++++++++++++++ script/Diagnostics/sig_process_wrapper.py | 31 +++ script/test/TestSigProcess.py | 30 +++ 5 files changed, 295 insertions(+), 10 deletions(-) create mode 100644 script/Diagnostics/sig_process.py create mode 100644 script/Diagnostics/sig_process_wrapper.py create mode 100644 script/test/TestSigProcess.py diff --git a/devices/CurrentCamera.properties b/devices/CurrentCamera.properties index 94dbd15..597ed2f 100644 --- a/devices/CurrentCamera.properties +++ b/devices/CurrentCamera.properties @@ -1,16 +1,16 @@ -#Wed May 03 16:50:12 CEST 2017 +#Thu May 04 11:39:44 CEST 2017 colormap=Flame colormapAutomatic=true colormapMax=800.0 colormapMin=0.0 -flipHorizontally=false +flipHorizontally=true flipVertically=false grayscale=false -imageHeight=1040 -imageWidth=1392 +imageHeight=2160 +imageWidth=2560 invert=false -regionStartX=0 -regionStartY=0 +regionStartX=1 +regionStartY=1 rescaleFactor=1.0 rescaleOffset=0.0 roiHeight=-1 @@ -21,9 +21,9 @@ rotation=0.0 rotationCrop=false scale=1.0 serverURL=localhost\:10000 -spatialCalOffsetX=-50.035945363048164 -spatialCalOffsetY=-50.04812319538017 -spatialCalScaleX=-1.0 -spatialCalScaleY=-1.0 +spatialCalOffsetX=-1279.0 +spatialCalOffsetY=-1079.0 +spatialCalScaleX=-53.020711215318485 +spatialCalScaleY=-53.02454840203798 spatialCalUnits=mm transpose=false diff --git a/script/Diagnostics/WireScan.py b/script/Diagnostics/WireScan.py index d19d77f..0aa84f7 100644 --- a/script/Diagnostics/WireScan.py +++ b/script/Diagnostics/WireScan.py @@ -164,6 +164,7 @@ finally: print "Closing stream" st.close() + # save the entry in the logbook if do_elog: if get_option("Generated data file:\n" + get_exec_pars().path +"\n\nSave to ELOG?", "YesNo") == "Yes": diff --git a/script/Diagnostics/sig_process.py b/script/Diagnostics/sig_process.py new file mode 100644 index 0000000..2234d88 --- /dev/null +++ b/script/Diagnostics/sig_process.py @@ -0,0 +1,223 @@ +# Signal processing functions +import copy +import numpy as np +import scipy.optimize + + +# Noise evaluation #### +def noise_evaluation(noise: np.ndarray): + """ + Returns offset and standard deviation of the noise signal. + :param noise: noise signal + :return: Tuple of (noise_offset, noise_sigma)) + """ + + return np.mean(noise), np.std(noise) + + +# Processing of BLM signal #### +def blm_remove_spikes(x): + """ + Returns new array with removed 1 and 2 point spikes. + :param x: input array + :return: output array + """ + + x = copy.copy(x) + if x.size > 5: # Must have enough sample points + d_x = x[1:] - x[:-1] + maximum = x.max() - x.min() + + for i in range(x.size-3): + if d_x[i+1] > 0.5 * maximum and d_x[i+2] < -0.5 * maximum: + # 1 point spike + x[i+2] = (x[i:i+5].sum() - x[i+2])/4 + + if d_x[i+1] > 0.5 * maximum and d_x[i+2] >= -0.5 * maximum: + # 2 point spikes + x[i+2] = (x[i:i+2].sum() + x[i+4])/3 + x[i+3] = (x[i+1] + x[i+4:i+6].sum())/3 + + # Handle edge points + if d_x[-1] > 0.5 * maximum: + x[-1] = x[-3:-1].sum() / 2 + + if d_x[-2] > 0.5 * maximum: + x[-2] = (x[-1] + x[-4:-2].sum()) / 3 + if d_x[0] < -0.5 * maximum: + x[0] = x[1:3].sum() / 2 + if d_x[1] < -0.5*maximum: + x[1] = (x[0] + x[2:4].sum()) / 3 + return x + + +def blm_normalize(x, q): + """ + Returns normalized BLM signal using charge from the BPM + :param x: BLM signal + :param q: BPM charge + :return: output array (same instance as input array) + """ + x = copy.copy(x) + return x / q + + +# Processing of position data #### +def motor_to_wire_cs(pos, type: str = 'u', center_pos: float = 0.0): + """ + Map wire position from motor coordinates (encoder readout) to wire coordinates. center_pos should be motor position + when wire is in the middle of the pipe. + + :param pos: position in motor coordinates + :param type: type of wire 'x', 'y' or 'u' (under 45 degrees to the pipe horizon) + :param center_pos: wire position in motor coordinates when in the middle of the pipe + :return: position in wire coordinates + """ + + # u: normal positions, like GARAGE, POSITION, or any wires under 45 degrees to pipe horizon + # x: which crosses pipe vertically (scans in x coordinate) + # y: which crosses pipe horizontally (scans in y coordinate) + w_scale = {'u': 1, 'x': -np.sqrt(2), 'y': np.sqrt(2)} + + return (pos-center_pos)/w_scale[type] + + +def remove_beam_jitter(pos: np.ndarray, bpm1: np.ndarray, bpm2: np.ndarray, d_b1_w=1, d_w_b2=1): + """ + This is a temporary solution which works only for first WSC on SwissFEL. + !!TODO: In future this method must implement correction using online model. + + :param pos: array of positions for 1 profile + :param bpm1: array of bpm in front of wsc readings (x or y) for 1 profile + :param bpm2: array of after wsc bpm readings (x or y) for 1 profile + :param d_b1_w: Distance between bpm1 and wsc + :param d_w_b2: Distance between wsc and bpm2 + :return: array of corrected positions + """ + + beam_position_at_wsc = (bpm1 * d_b1_w + bpm2 * d_w_b2) / (d_b1_w + d_w_b2) + + # Reference centroid position is mean value of all measured values + beam_jitter = beam_position_at_wsc - beam_position_at_wsc.mean(0) + + # Correct wire position for beam jitter + return pos + beam_jitter + + +# Profile statistics #### +def profile_gauss_stats(x, y, off=None, amp=None, com=None, sigma=None): + if off is None: + off = y.min() # good enough starting point for offset + + if com is None: + com = x[y.argmax()] + + if amp is None: + amp = y.max() - off + + # For normalised gauss curve sigma=1/(amp*sqrt(2*pi)) + if sigma is None: + surface = np.trapz((y-off), x=x) + sigma = surface / (amp * np.sqrt(2 * np.pi)) + try: + popt, pcov = scipy.optimize.curve_fit(gauss_fn, x, y, p0=[off, amp, com, sigma], method='lm') + popt[3] = abs(popt[3]) # sigma should be returned as positive + except Exception as e: + print("Gauss fitting not successful.\n" + str(e)) + popt = [off, amp, com, abs(sigma)] + + return popt + + +def gauss_fn(x, a, b, c, d): + return a + b * np.exp(-(np.power((x - c), 2) / (2 * np.power(d, 2)))) + + +def profile_rms_stats(x, y, noise_std=0, n_sigma=3.5): + """ + Does the center of mass and RMS calculation on given profile. + + :param x: data of x axis (wire position) + :param y: data of y axis (blm current) + :param noise_std: standard deviation of noise + :param n_sigma: multiplication of window size (n_sigma*sigma) where sigma is first estimation of sigma + :return: rms of signal + """ + + if x.shape == y.shape: + # Windowing the signal + # 1. All points bellow +sigma are treated as noise and ignored + # 2. If point is higher than +sigma and neighbour points are noise, then it is also treated as noise and ignored + # 3. All sets of multiple non noise points in a row are treated as signal_candidates + signal_flag = False + signal_start = 0 + signal_candidates = list() + for i in range(y.size): + if y[i] < noise_std: + # Not part of the signal. Check if this is the end of signal candidate. + if signal_flag and i-signal_start > 1: + signal_candidates.append((signal_start, i-1)) + signal_flag = False + + elif not signal_flag: + # Start of potential new signal candidate + signal_flag = True + signal_start = i + + # Take care of last signal candidate if last point is higher than noise (signal_flag is still True) + if signal_flag and (y.size - signal_start - 1) > 1: + signal_candidates.append((signal_start, y.size-1)) + + if len(signal_candidates) == 0: + # Only noise, no signal. Cannot calculate rms + print("RMS calculation not successful. No signal detected.\n") + return None, None + + # Evaluate candidates + candidates_val = [y[s[0]:s[1]].sum() for s in signal_candidates] + signal_boundaries = signal_candidates[candidates_val.index(max(candidates_val))] + + # Calculate RMS values (at this point they are underestimated, but sigma is then used to extend the signal + # region of interest which gives a better result (biggest surface) + com_estimate, rms_estimate = com_rms(x[signal_boundaries[0]:signal_boundaries[1]], + y[signal_boundaries[0]:signal_boundaries[1]]) + + # Indexes of useful signal (-1 is here because original algorithm calculates this a bit different and + # because it does the floor() of calculated function + + # Define direction of scan + if x[0] < x[-1]: # increasing position + low_boundary = np.argmax(x >= com_estimate - n_sigma * rms_estimate)-1 + high_boundary = np.argmax(x >= com_estimate + n_sigma * rms_estimate)-1 + + else: # decreasing position + high_boundary = np.argmax(x <= com_estimate - n_sigma * rms_estimate)-1 + low_boundary = np.argmax(x <= com_estimate + n_sigma * rms_estimate)-1 + + if high_boundary <= 0: + high_boundary = x.size + + if low_boundary < 0: + low_boundary = 0 + + if (high_boundary-low_boundary) < 4: + print("RMS calculation not successful. Resolution to low.\n") + return None, None + + # Calculate final RMS and centroid on ROI data + if high_boundary == x.size: + com, rms = com_rms(x[low_boundary:], y[low_boundary:]) + else: + com, rms = com_rms(x[low_boundary:high_boundary], y[low_boundary:high_boundary]) + + return com, rms + else: + return None, None + + +def com_rms(x, y): + # Centre of mass and rms + com = (x * y).sum() / y.sum() + com2 = (x * x * y).sum() / y.sum() + rms = np.sqrt(np.abs(com2 - com * com)) + return com, rms \ No newline at end of file diff --git a/script/Diagnostics/sig_process_wrapper.py b/script/Diagnostics/sig_process_wrapper.py new file mode 100644 index 0000000..62b9020 --- /dev/null +++ b/script/Diagnostics/sig_process_wrapper.py @@ -0,0 +1,31 @@ +from jeputils import * + +MODULE = "Diagnostics/sig_process" + +def noise_evaluation(noise): + return call_jep(MODULE, "noise_evaluation", [to_npa(noise),]) + +def blm_remove_spikes(x): + ret = call_jep(MODULE, "blm_remove_spikes", [to_npa(x),]) + return None if ret is None else ret.data + +def blm_normalize(x, q): + ret = call_jep(MODULE, "blm_normalize", [to_npa(x), q]) + return None if ret is None else ret.data + +def motor_to_wire_cs(pos, ctype = 'u', center_pos = 0.0): + return call_jep(MODULE, "motor_to_wire_cs", [pos, ctype, center_pos ]) + +def remove_beam_jitter(pos, bpm1, bpm2, d_b1_w=1, d_w_b2=1): + ret = call_jep(MODULE, "remove_beam_jitter", [to_npa(pos),to_npa(bpm1), to_npa(bpm2), d_b1_w, d_w_b2 ]) + return None if ret is None else ret.data + +def profile_gauss_stats(x, y, off=None, amp=None, com=None, sigma=None): + ret = call_jep(MODULE, "profile_gauss_stats", [to_npa(x), to_npa(y), off, amp, com, sigma]) + return None if ret is None else ret.data + +def profile_rms_stats(x, y, noise_std=0, n_sigma=3.5): + return call_jep(MODULE, "profile_rms_stats", [to_npa(x), to_npa(y), noise_std, n_sigma]) + + + \ No newline at end of file diff --git a/script/test/TestSigProcess.py b/script/test/TestSigProcess.py new file mode 100644 index 0000000..dbe3830 --- /dev/null +++ b/script/test/TestSigProcess.py @@ -0,0 +1,30 @@ +run("Diagnostics/sig_process_wrapper") + +#Loading test data +#p=get_plots("Wire Scan")[0] +#s=p.getSeries(0) +#indexes = sorted(range(len(s.x)),key=lambda x:s.x[x]) +#x,y = [s.x[x] for x in indexes], [s.y[x] for x in indexes] +x=[-200.30429237268825, -200.2650700434188, -200.22115208318002, -199.9457671375377, -199.86345548879072, -199.85213073174933, -199.35687977133284, -199.13811861090275, -197.97304970346386, -197.2952215624348, -195.09076092936948, -192.92276048970703, -191.96871876227698, -189.49577852322938, -187.9652790409825, -183.63756456925222, -180.04899765472996, -178.43839623242422, -174.07311671294445, -172.0410133577918, -165.90824309893102, -160.99771795989466, -159.30176653939253, -154.27688897558514, -152.0854103810786, -145.75652847587313, -140.80843828908465, -139.23982133191495, -134.27073891256106, -132.12649284133064, -125.95947209775511, -121.00309550337462, -119.26736932643232, -114.2706655484383, -112.07393889578914, -105.72295990367157, -100.8088439880125, -99.2034906238494, -94.30042325164636, -92.15010048151461, -85.92203653534293, -81.03913275494665, -79.27412793784428, -74.33487658582118, -72.06274362408762, -65.76562628131825, -60.91255356825276, -59.20334389560392, -54.33286972659312, -52.19387171350535, -45.94978737932291, -41.03014719193582, -39.301602568238906, -34.35572209014114, -32.04464301272608, -25.8221033382824, -20.922074315528747, -19.21590299233186, -14.31090212502093, -12.217203140101386, -5.9283722049240435, -0.9863587170369246, 0.7408048387279834, 5.71126832601389, 7.972628957879352, 14.204559894256546, 19.11839959633025, 20.8218087836657, 25.678748486941828, 27.822718344586864, 34.062659474970715, 38.9745656819391, 40.77409719734158, 45.72080631619803, 47.974156754056835, 54.23453768983539, 59.12020360609568, 60.77306570712026, 65.70734521458867, 67.8344660434617, 74.03187028154134, 78.96532114824849, 80.76070945985495, 85.74802197591286, 87.9140889204674, 94.18082276873524, 99.25790470037091, 100.68454787413205, 105.7213026221542, 107.79483801526698, 113.99555681638138, 119.0707052529143, 120.72715813056156, 125.77551384921307, 127.91257836719551, 134.2011330887875, 139.23043006997628, 140.71673537840158, 145.76288138835983, 147.80216629676042, 154.06420451405637, 159.0846626604798, 160.76183155710717, 165.73699067536242, 167.9265357747636, 173.96705069576544, 178.2522282751915, 179.9042617354548, 183.54586165856657, 185.23269803071796, 189.41678143751972, 191.87149157986588, 192.8741468985015, 195.0241934550453, 195.966634211846, 197.9821647518146, 198.99006812859284, 199.33202054855676, 199.91897441965887, 200.11536227958896, 200.22280936469997, 200.25181179127208] +y=[11.0, 6.0, 8.0, 5.0, 11.0, 7.0, 18.0, 11.0, 12.0, 10.0, 8.0, 6.0, 16.0, 4.0, 12.0, 9.0, 15.0, 14.0, 8.0, 20.0, 15.0, 8.0, 9.0, 11.0, 13.0, 12.0, 13.0, 15.0, 13.0, 20.0, 10.0, 7.0, 17.0, 11.0, 20.0, 13.0, 13.0, 23.0, 14.0, 10.0, 17.0, 15.0, 20.0, 16.0, 14.0, 13.0, 18.0, 22.0, 9.0, 20.0, 12.0, 14.0, 17.0, 19.0, 14.0, 14.0, 23.0, 19.0, 15.0, 20.0, 20.0, 21.0, 20.0, 23.0, 22.0, 15.0, 10.0, 17.0, 21.0, 15.0, 23.0, 23.0, 25.0, 18.0, 16.0, 21.0, 22.0, 16.0, 16.0, 14.0, 19.0, 20.0, 18.0, 20.0, 23.0, 13.0, 16.0, 20.0, 25.0, 15.0, 15.0, 17.0, 22.0, 26.0, 19.0, 30.0, 25.0, 17.0, 17.0, 23.0, 16.0, 27.0, 21.0, 21.0, 26.0, 27.0, 21.0, 17.0, 20.0, 20.0, 21.0, 19.0, 25.0, 19.0, 13.0, 23.0, 20.0, 20.0, 18.0, 20.0, 19.0, 25.0] +y[30] = 100 + +#Calling wrapper methods +rem= blm_remove_spikes(y) +print "Noise: ", noise_evaluation(rem) +print "Wire X: ", motor_to_wire_cs(20000, 'x', 2000) +nor = blm_normalize(rem, 2.0) +s2 = [i + random.random()*10 for i in nor] +jit = remove_beam_jitter(x, nor, s2, d_b1_w=1, d_w_b2=1) +[com, rms] = profile_rms_stats(x, nor,noise_std=0, n_sigma=3.5) +print "Rms: ", [com, rms] +[off, amp, com, sigma] = profile_gauss_stats(x, nor, off=None, amp=None, com=None, sigma=None) +print "Gauss: ", [off, amp, com, sigma] +from mathutils import Gaussian +f = Gaussian(amp, com, sigma) +gauss = [f.value(i)+off for i in x] + +#Plotting results +plot([y,rem, nor, gauss, s2, jit], ["signal", "rem", "nor", "gauss", "s2", "jit"], xdata = x, title="Fit") + +