AareDAQ/scripts/knife_edge_scan.ipynb

Import Packages

In [1]:

from aaredaq.devices import BeamlineDevices
from aaredaqlib.beamline import MXBeamline

In [2]:

from epics import PV
import numpy as np
import time
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
from scipy.signal import peak_widths, find_peaks

Config

In [3]:

d = BeamlineDevices(MXBeamline.X06DA)
b = MXBeamline.X06DA
beamline = b.value.upper()

Connecting TELL p-shell service at http://x06da-tell.psi.ch:22222 ...

Functions for fitting

In [4]:

def gaus(x, *p):
    A, mu, sigma = p
    return A * np.exp(-(x-mu)**2/(2.*sigma**2))

def sigmoid(x, *p):
    L, x0, k, c = p
    return L / (1 + np.exp(-k * (x -x0))) + c

def unpack_data(data):
    x_vals, y_vals = zip(*data)
    x_vals = list(x_vals)
    y_vals = list(y_vals)
    return x_vals, y_vals

def sigmoid_fit(x_vals, y_vals):
    L_guess = np.max(y_vals)
    b_guess = np.min(y_vals)
    x0_guess = x_vals[np.argmin(np.abs(y_vals - (L_guess + b_guess) / 2))]
    k_guess = 1 / (x_vals[-1] - x_vals[0])

    s0 = [L_guess, x0_guess, k_guess, b_guess]
    parameters, params_covariance = curve_fit(sigmoid, x_vals, y_vals, p0=s0)
    return parameters

def gauss_fit(x_vals, dydx):
    mean = x_vals[np.argmax(dydx)]
    sigma = (np.max(x_vals) - np.min(x_vals))/4
    A= np.max(dydx)
    dydx[np.isnan(dydx)] = 0

    p0 = [A, mean, sigma]
    parameters, covariance = curve_fit( gaus, x_vals, dydx, p0=p0)
    return parameters

def calc_fwhm(gauss_y, x_vals):
    peaks = find_peaks( gauss_y )
    fwhm_peak = peak_widths( gauss_y, peaks[0], rel_height=0.5 )

    fwhm_str = x_vals[int( round( fwhm_peak[2][0], 0 ))]
    fwhm_end = x_vals[int( round( fwhm_peak[3][0], 0 ))]

    fwhm = fwhm_str - fwhm_end
    fwhm_height = gauss_y[int( round( fwhm_peak[2][0], 0 ))]

    # find 1/e2 peak width
    full_peak = peak_widths( gauss_y, peaks[0], rel_height=0.865 )
    full_str = x_vals[int( round( full_peak[2][0], 0 ))]
    full_end = x_vals[int( round( full_peak[3][0], 0 ))]

    full = full_str - full_end
    full_height = gauss_y[int( round( full_peak[2][0], 0 ))]
    return fwhm, full

Function for performing knife_edge then fitting and plotting

In [41]:

def knife_edge(motor: str, dev: BeamlineDevices, beamline: str, steps: int = 40, start: float = 0.0, step_size: float = 0.002, offset: float = 0.0, sleep = 0.5):
    diode = PV(f"{beamline}-ES-BS:READOUT")#X06DA-ES-BS:READOUT
    dev.shutter = True
    time.sleep(2.0)

    data = []

    for i in range(steps): #40
        if motor == 'X':
            dev.gmx.move(start-offset+step_size*i, wait=True)
            time.sleep(sleep)
            print(f"{dev.gmx.value:.3f} {diode.get()}")
            data.append([dev.gmx.value, diode.get()])
        else:
            dev.gmy.move(start-offset+step_size*i, wait=True)
            time.sleep(sleep)
            print(f"{dev.gmy.value:.3f} {diode.get()}")
            data.append([dev.gmy.value, diode.get()])

    if motor == 'X':
        dev.gmx.move(start, wait=True)
    #else:
    #    dev.gmy.move(start, wait=True)

    dev.shutter = False
    return data

In [42]:

def fit_knife_edge(data, motor: str):

    x_vals, y_vals = unpack_data(data)

    plt.plot(x_vals, y_vals, label = 'edge_scan')
    plt.ylabel('Intensity (counts)')
    plt.xlabel('Motor position (mm)')

    params=sigmoid_fit(x_vals, y_vals)

    plt.plot(x_vals, y_vals)
    plt.plot(x_vals, sigmoid(x_vals, *params))
    plt.show()

    y_fit = sigmoid(x_vals, *params)

    dydx = np.gradient(y_fit, x_vals)

    if motor == 'X':
        dydx=-dydx

    parameters = gauss_fit(x_vals, dydx)

    x0_op = parameters[1]
    sigma_op = parameters[2]
    print(parameters)

    gauss_y = gaus(x_vals,*parameters)
    fwhm = np.abs(2*np.sqrt(2*np.log(2))*sigma_op)
    print(f"Manually calculated FWHM to {fwhm} mm")
    plt.plot(x_vals, dydx, label = 'derivative')
    plt.plot(x_vals, gauss_y,label='Gaussian fit',color ='orange')
    plt.fill_between(x_vals,gauss_y,color='orange',alpha=0.5)

    try:
        fwhm, full = calc_fwhm(gauss_y, x_vals)
        print( "Scipy calculated FWHM = {0} mm".format( fwhm ) )
        print( "Scipy calculated 1/e2 = {0} mm".format( full ) )
    except:
        print('scipy peak finding failed')

    plt.axvspan(x0_op+fwhm/2,x0_op-fwhm/2, color='green', alpha=0.75, lw=0, label='FWHM = {0} mm'.format(fwhm))
    plt.legend()
    plt.show()
    #plt.savefig(output_name+filename.split('/')[-1]+'FWHM_{0}.png'.format(fwhm))

Knife Edge scan in X

In [63]:

steps_x = 50
start_pos_x = -17.94
step_size_x = 0.002
offset_x = 0.0 # offset subtracted from start position
sleep = 1 # wait after move before recording diode
x_data=knife_edge(motor='X', dev=d, beamline=beamline, steps=steps_x, start=start_pos_x, step_size=step_size_x, offset=offset_x, sleep = sleep)
fit_knife_edge(x_data, motor='X')

-17.940 8.959981e-06
-17.938 8.944953e-06
-17.936 8.937235e-06
-17.934 8.946256e-06
-17.932 8.870569e-06
-17.930 8.943733e-06
-17.928 8.932527e-06
-17.926 8.931002e-06
-17.924 8.912179e-06
-17.922 8.929313e-06
-17.920 8.917321e-06
-17.918 8.937232e-06
-17.916 8.892818e-06
-17.914 8.824038e-06
-17.912 8.808735e-06
-17.910 8.837087e-06
-17.908 8.671096e-06
-17.906 8.606978e-06
-17.904 8.432655e-06
-17.902 8.189082e-06
-17.900 8.250342e-06
-17.898 7.650508e-06
-17.896 7.203671e-06
-17.894 6.735723e-06
-17.892 6.1932e-06
-17.890 6.269931e-06
-17.888 5.137038e-06
-17.886 4.552743e-06
-17.884 3.875157e-06
-17.882 3.34356e-06
-17.880 3.422219e-06
-17.878 2.437191e-06
-17.876 2.017403e-06
-17.874 1.662639e-06
-17.872 1.323577e-06
-17.870 1.346216e-06
-17.868 9.006973e-07
-17.866 7.306933e-07
-17.864 5.915996e-07
-17.862 4.870485e-07
-17.860 4.98184e-07
-17.858 3.496483e-07
-17.856 2.945638e-07
-17.854 2.461977e-07
-17.852 2.114641e-07
-17.850 2.139551e-07
-17.848 1.598115e-07
-17.846 1.397988e-07
-17.844 1.199971e-07
-17.842 1.045175e-07

[ 2.92142522e-04 -1.78855968e+01  1.17269351e-02]
Manually calculated FWHM to 0.02761482178449844 mm
Scipy calculated FWHM = -0.02800645446777139 mm
Scipy calculated 1/e2 = -0.04800543212890318 mm

In [ ]:

fit_knife_edge(x_data, motor='X')

Knife edge scan in Y

In [68]:

steps_y = 40
start_pos_y = -0.20
step_size_y = 0.002
offset_y = 0.00
sleep = 1.0 # wait after move before recording diode
y_data=knife_edge(motor='Y', dev=d, beamline=beamline, steps=steps_y, start=start_pos_y, step_size=step_size_y, offset=offset_y, sleep=sleep)
fit_knife_edge(y_data, motor='Y')

-0.200 2.360645e-08
-0.198 2.811506e-08
-0.196 2.810113e-08
-0.194 3.39931e-08
-0.192 4.061688e-08
-0.190 6.625245e-08
-0.188 8.658314e-08
-0.186 1.106975e-07
-0.184 1.428539e-07
-0.182 1.863939e-07
-0.180 1.863939e-07
-0.178 2.407869e-07
-0.176 2.862987e-07
-0.174 3.589597e-07
-0.172 4.629766e-07
-0.170 7.34634e-07
-0.168 9.591903e-07
-0.166 1.226978e-06
-0.164 1.638941e-06
-0.162 2.070232e-06
-0.160 2.070232e-06
-0.158 2.719557e-06
-0.156 3.40311e-06
-0.154 4.3336e-06
-0.152 5.30817e-06
-0.150 7.05408e-06
-0.148 7.882228e-06
-0.146 8.391648e-06
-0.144 8.645698e-06
-0.142 8.800723e-06
-0.140 8.833474e-06
-0.138 8.667967e-06
-0.136 8.724034e-06
-0.134 8.751249e-06
-0.132 8.949502e-06
-0.130 8.983748e-06
-0.128 8.932714e-06
-0.126 8.97258e-06
-0.124 8.960077e-06
-0.122 9.000123e-06

[ 0.00048058 -0.15422773  0.00713343]
Manually calculated FWHM to 0.016797947350266428 mm
Scipy calculated FWHM = -0.015900000000000247 mm
Scipy calculated 1/e2 = -0.02798400000000001 mm

In [ ]:

fit_knife_edge(y_data, motor='Y')