Files
pxii_bec/pxii_bec/macros/APscripts.py
T
x10sa 7147bc1c84
CI for pxii_bec / test (push) Successful in 36s
Add beamline states
2026-04-21 10:46:29 +02:00

1507 lines
43 KiB
Python
Executable File

################################
### Scripts/"Macros" for BEC:
################################
# Python scripts for running scan plans in BEC are to be stored in
# ~/bec/scripts/
import numpy as np
import math
import matplotlib.pyplot as plt
################################
### BASIC BL PROGRAMS
#####################
### a2e
#####################
def a2e(a, *hkl):
import math
# possible keywords: # dictionary in python
# DEG = deg
# LN = ln
i = 0
count = 0
h = [1, 1, 1]
iln = 0
ideg = 0
while i < 3 and count < len(hkl):
for var in hkl:
count = count + 1
# print(var)
if type(var) == int:
print("var int ", var, i)
h[i] = var
i = i + 1
if "ln" in hkl:
iln = 1
if "deg" in hkl:
ideg = 1
d0 = 2 * 5.43102 * (1 - 2.4e-4 * iln) / np.sqrt(h[0] ** 2.0 + h[1] ** 2.0 + h[2] ** 2.0)
if ideg or (a > 1):
a = math.radians(a) # *math.pi/180.
en = 12.39842 / d0 / math.sin(a)
return en
#####################
### angle
#####################
def angle(e, *hkl):
import math
# possible keywords: # dictionary in python
# DEG = deg
# LN = ln
i = 0
count = 0
h = [1, 1, 1]
iln = 0
ideg = 0
while i < 3 and count < len(hkl):
for var in hkl:
count = count + 1
# print(var)
if type(var) == int:
print("var int ", var, i)
h[i] = var
i = i + 1
if "ln" in hkl:
iln = 1
if "deg" in hkl:
ideg = 1
d0 = 2 * 5.43102 * (1 - 2.4e-4 * iln) / np.sqrt(h[0] ** 2.0 + h[1] ** 2.0 + h[2] ** 2.0)
a = math.asin(12.39842 / d0 / e)
if ideg:
a = math.degrees(a) # *180./math.pi
return a
#####################
### setenergy
#####################
#####################
def sete(e, rock=1):
bragg_ang = angle(e, "ln") * 1000.0
umv(dev.dcm_bragg, bragg_ang)
setu19(e, detune=1)
if rock:
rock()
return
#####################
### getenergy
#####################
def getenergy():
ang_bragg_rb = dev.dcm_bragg.read()["dcm_bragg"]["value"]
# dev.dcm_bragg.user_readback.get()
ang_bragg = ang_bragg_rb / 1000.0 # mrad to rad
energy = a2e(ang_bragg, "ln")
return round(energy, 5)
#####################
### ROCK
#####################
def rock(**kwargs):
from lmfit.models import LinearModel, GaussianModel # LorentzianModel
import matplotlib.pyplot as plt
dock_area = bec.gui.new()
wr = dock_area.new().new(bec.gui.available_widgets.Waveform)
# width of rocking curve of perfect xtal: at 20 keV: 14 urad == 0.0008 deg
# keywords:
# silent (no terminal output)
# ffail (success? repeat ?)
# time
# dth - change range of rocking curve (e.g.,if far off)
### axis [0/1] 0: Bragg ; 1: pitch
if "axis" in list(kwargs):
ax = kwargs["axis"]
print("Rocked MonoAxis is dcm_bragg")
else:
ax = 0
# print('Rocked MonoAxis is Th1')
print("Rocked MonoAxis is dcm_pitch")
### waiting time
if "time" in list(kwargs):
time = kwargs["time"]
else:
time = 0.3
print("Time is default of 0.3 s")
print("Time is %s s" % time)
e = getenergy()
### width of rocking curve
if "dtheta" in list(kwargs):
dx = kwargs["dtheta"]
else:
# dx = 0.05/e ### width: 0.1/6 to 0.1/30 ~~ 0.017 to 0.0033 ###### in deg
dx = 0.1 # in mrad## width: 0.1/6 to 0.1/30 ~~ 0.017 to 0.0033
dx = 0.6 / e * 2
print("Rocking width is %s mrad" % dx)
det = dev.lu_bpmsum
mot = dev.dcm_pitch
motb = dev.dcm_bragg
mot_rb = mot.user_readback.get()
motb_rb = motb.user_readback.get()
# or pitch0 = dev.dcm_pitch.read()['dcm_pitch']['value']
pos0 = mot_rb
print("Position of bragg angle is %s mrad" % motb_rb)
print("Position of bragg pitch is %s mrad" % mot_rb)
if ax == 1:
mot = motb
pos0 = motb_rb
s = scans.line_scan(mot, -dx, dx, steps=50, exp_time=time, relative=True)
# md = scan.metadata["bec"]
wr.title = f"RockingScan at Energy of {e}"
wr.plot(x_name=mot.name, y_name=det.name) ##set names/axes first !
wr.x_label = mot.name
wr.y_label = det.name
if ax == 0:
data_x = s.scan.live_data.dcm_pitch.dcm_pitch.val
else:
data_x = s.scan.live_data.dcm_bragg.dcm_bragg.val
data_y = s.scan.live_data.lu_bpmsum.lu_bpmsum.val
wr.plot(data_x, data_y)
# fit Gaussian + lin BG
peak = GaussianModel()
background = LinearModel()
model = peak + background
maxy = max(data_y)
indmax = np.argmax(data_y)
xm = data_x[indmax]
p = model.make_params(amplitude=maxy, center=xm)
g = model.fit(data_y, p, x=data_x)
print(f'Center: {g.params["center"].value:.5f} FWHM: {g.params["fwhm"].value:.5f}')
plt.ion()
plt.figure(1)
plt.plot(data_x, data_y, "*")
plt.plot(data_x, g.best_fit)
plt.show()
res_center = g.params["center"].value
diff_movement = abs(res_center - pos0)
if diff_movement < 0.25:
umv(mot, res_center)
else:
print("Fit failed, no movement")
return
################################
### FITTING PROGRAMS
#############################
### fit and plot of history
#############################
def justfit(data_x, data_y, model="gauss", ibg=0):
"""
just fit a fct of Gauss, Voigt or Lorentzian type,
Examples:
gfit, xmax = justfit(data_x, data_y, model = "gauss", ibg =1) : Gauss, lin BG
gfit, xmax = justfit(data_x, data_y, model = "voigt", ibg =2) : Voigt, quadrat BG
gfit, xmax = justfit(data_x, data_y, model = "lorentz", ibg =0) : Lorentzian, no BG
"""
from lmfit.models import LinearModel, GaussianModel, VoigtModel, QuadraticModel, LorentzianModel
import matplotlib.pyplot as plt
peak = GaussianModel()
if model == "voigt":
peak = VoigtModel()
if model == "lorentz":
peak = LorentzianModel()
model = peak
if ibg:
# add a BG
background = LinearModel()
if ibg > 1:
background = QuadraticModel()
model = peak + background
### define data_x and data_y:
maxy = max(data_y)
indmax = np.argmax(data_y)
xm = data_x[indmax]
sigma = 1.0
gamma = 0.2 # blurring/widening of the sigma ; the larger, the more of a Lorentzian profile
## possible, but not necessarily better ... :
# sl = 0.00001 * (data_y[-2]-data_y[2])/(data_x[-2]-data_x[2]) # small slope
# ic = min(data_y)
# print(f'slope and intercept initial conds= {sl, ic}')
##
print("maxy, indmax, xm = ", maxy, indmax, xm)
p = model.make_params(amplitude=maxy, center=xm)
# other possibilities ... :
# p = model.make_params(amplitude = maxy, center=xm , slope=0, intercept = ic)
# p = model.make_params(amplitude = max(data_y), center=xm, slope=0, intercept = min(data_y))
p["center"].set(min=min(data_x), max=max(data_x))
p["sigma"].set(min=0, max=(max(data_x) - min(data_x)) / 2.0)
# other possibilities ... :
# p['amplitude'].set(min=min(data_y), max=max(data_y))
# p['intercept'].set(max=(max(data_y)-min(data_y))/2.)
g = model.fit(data_y, p, x=data_x)
# diagnostics
# print(f'Gfit: {g.params}')
print(f'Center of {model} fit: {g.params["center"].value:.5f}')
print(f'Sigma: {g.params["sigma"].value:.5f}')
print(f'FWHM: {g.params["fwhm"].value:.5f}')
print(f"Position of maximum: {xm}")
return g, xm
def fit_plothist(hindex: int, signal_name: str, model="gauss", ibg=0):
"""
Get data for a completed scan from the BEC history, (gaussian, voigt, lorentzian + zero, lin, quadrat BG) fit,
and plot all the results.
Args:
history_index (int): scan to fetch, e.g. -1 for the most recent scan
signal_name (str): the monitor
Example:
gcen, xm , data_x, data_y = fit_plothist(-10, "lu_bpmsum",model="gauss", ibg=2)
"""
from lmfit.models import LinearModel, GaussianModel, VoigtModel, QuadraticModel, LorentzianModel
import matplotlib.pyplot as plt
h = bec.history[hindex]
md = h.metadata["bec"]
scanvar = list(md["args"].keys())[0] # string, returns the variable of the last performed scan
# data = h.devices[device_name][signal].read()["value"]
# data_x = h.devices.dcm_pitch.dcm_pitch.read()["value"]
# data_y = h.devices.lu_bpmsum.lu_bpmsum.read()["value"]
# data_x = h.devices[device_name][device_name].read()["value"]
data_xd = md["positions"] # last scan knows which device ... however, double array
data_x = data_xd.flatten()
print("dx = ", data_x)
data_y = h.devices[signal_name][signal_name].read()["value"]
print("dy = ", data_y)
g, xm = justfit(data_x, data_y, model=model, ibg=ibg)
plt.ion()
plt.figure()
plt.plot(data_x, data_y, "*")
plt.plot(data_x, g.best_fit)
plt.xlabel(scanvar)
plt.ylabel(signal_name)
plt.show()
gcen = g.params["center"].value
return gcen, xm, data_x, data_y
#########################
### fit and plot ONLY
#########################
def fit_plot(data_x, data_y, model="gauss", ibg=1, fitrange=0, fitclick=0):
"""
Get data for a completed scan, gaussian, voigt or lorentzian fit,
default = gauss, plus linear background
and plot all the results;
possibly select the range before [fitrange = 1], you will be asked to enter the x-range
OR choose interactively [[fitclick = 1], you will be asked to click
Args:
datax , datay
Example:
fit_plot(dx,dy, model='voigt', fitrange=0, fitclick=1)
"""
from lmfit.models import LinearModel, GaussianModel, VoigtModel, LorentzianModel
import matplotlib.pyplot as plt
import numpy as np
print("dy = ", data_y)
####### check the order of the data: lower --> upper ? if not, reverse all the indices!
sz = len(data_x)
if data_x[sz - 1] < data_x[0]:
# reverse indices
data_x = data_x[::-1]
data_y = data_y[::-1]
if fitrange:
x1, x2 = input("Please enter the x-range for the fit as x1, x2 ").split(",")
x1 = float(x1)
x2 = float(x2)
closest_low = np.abs(np.array(data_x) - x1).argmin()
closest_high = np.abs(np.array(data_x) - x2).argmin()
print("low_ind", closest_low)
print("high_ind", closest_high)
data_x = data_x[closest_low:closest_high]
data_y = data_y[closest_low:closest_high]
print("dx=", data_x)
print("dy=", data_y)
#####
peak = GaussianModel()
if model == "voigt":
peak = VoigtModel()
if model == "lorentz":
peak = LorentzianModel()
if ibg:
background = LinearModel()
model = peak + background
### define data_x and data_y:
maxy = max(data_y)
indmax = np.argmax(data_y)
xm = data_x[indmax]
sigma = 1.0
gamma = 0.2 # blurring/widening of the sigma ; the larger, the more of a Lorentzian profile
print("maxy, indmax, xm = ", maxy, indmax, xm)
# p = model.make_params(
# amplitude=max(data_y), center=xm, slope=0, intercept=min(data_y)
# )
p = model.make_params(amplitude=maxy, center=xm)
p["center"].set(min=min(data_x), max=max(data_x))
p["sigma"].set(min=0, max=(max(data_x) - min(data_x)) / 2.0)
g = model.fit(data_y, p, x=data_x)
# diagnostics
# print(f'Gfit: {g.params}')
print(f'Center of {model} fit: {g.params["center"].value:.5f}')
print(f'Sigma: {g.params["sigma"].value:.5f}')
print(f'FWHM: {g.params["fwhm"].value:.5f}')
print(f"Position of maximum: {xm}")
#####
plt.ion()
plt.figure()
plt.plot(data_x, data_y, "*")
plt.plot(data_x, g.best_fit)
plt.xlabel("x-axis")
plt.ylabel("signal")
plt.show()
if fitclick:
print("Please define the x-range of the fit via click of graph")
print("First click for start: ")
x1 = plt.ginput(1)
print("Second click for end: ")
x2 = plt.ginput(1)
print("x1 = ", x1)
print("x2 = ", x2)
x1 = float(x1[0][0]) ## first element of the tupel in the list
x2 = float(x2[0][0])
print("x1,x2 = ", x1, x2)
closest_low = np.abs(data_x - x1).argmin()
closest_high = np.abs(data_x - x2).argmin()
print("low_ind", closest_low)
print("high_ind", closest_high)
data_x = data_x[closest_low:closest_high]
data_y = data_y[closest_low:closest_high]
print("dx=", data_x)
print("dy=", data_y)
g, xm = justfit(data_x, data_y, model=model, ibg=ibg)
plt.figure()
plt.plot(data_x, data_y, "*")
plt.plot(data_x, g.best_fit)
plt.xlabel("x-axis")
plt.ylabel("signal")
plt.show()
gcen = g.params["center"].value
return gcen, xm
####################################
### just save data
####################################
def save_data(hindex: int, device_name: str, signal_name: str):
"""
Get data for a completed scan from the BEC history,
and plot all the results, store in a CSV file
Args:
history_index (int): scan to fetch, e.g. -1 for the most recent scan
device_name (str): the device for which to get the monitoring data
signal_name (str): the monitor
"""
h = bec.history[hindex]
md = h.metadata["bec"]
print("Scan argument was ", md["args"])
data_x = h.devices[device_name][device_name].read()["value"]
data_y = h.devices[signal_name][signal_name].read()["value"]
plt.ion()
plt.figure()
plt.plot(data_x, data_y, ".")
plt.title(f"Scan of {device_name}")
plt.xlabel(f"{device_name}")
plt.ylabel(f"{signal_name} / AU")
plt.show()
ans = "n"
ans = input("Store data in csv file? y/n ")
if ans == "y":
dirname = "/sls/x10sa/config/commissioning/Data/"
# writing output to simple data file for later analysis:
combined = np.column_stack((data_x, data_y))
filename = dirname + "Scan" + str(hindex) + device_name + ".txt"
with open(filename, "w") as f:
np.savetxt(f, combined, delimiter=",", fmt="%5f")
return data_x, data_y
#########################################
### just retrieve the saved data from csv
########################################
def read_data(filename: str):
"""
Get data stored in a CSV file
Args:
filename (str): the csv file, eg, of a Scan
e.g., dat = read_data("/home/e18747/SLS2/Data/gaps/gaps10.txt")
"""
# dirname = '/home/gac-x10sa/Data/'
# read output to simple data file for later analysis:
with open(filename, "r") as f:
combined_data = np.loadtxt(f, delimiter=",") # no header and fmt for load
rows, cols = combined_data.shape
print("No of Rows =", rows)
print("No of Colums =", cols)
ind1 = 0
ind2 = 1
if cols > 2:
print(
"only first 2 colums (data[:, 0] and data[:, 1]) are plotted, please consider the other columns as well!"
)
index1 = 0
index2 = 1
index1, index2 = input(
"please enter the columns to be plotted vs first one [eg: 1,2] : "
).split(",")
ind1 = int(index1)
ind2 = int(index2)
if ind2 > cols - 1:
print(" no valid index")
data_x = combined_data[:, ind1]
data_y = combined_data[:, ind2]
plt.ion()
plt.figure()
plt.plot(data_x, data_y, ".")
plt.title(f"Scan from {filename}")
plt.xlabel(f"data column {ind1}")
plt.ylabel(f"data column {ind2}")
plt.show()
return combined_data
####################################
### save and plot the long gapscans
####################################
def save_plot_gaps(hindex: int, device_name: str, signal_name: str):
"""
Get data for a completed scan from the BEC history,
and plot all the results, store in a CSV file
Args:
history_index (int): scan to fetch, e.g. -1 for the most recent scan
device_name (str): the device for which to get the monitoring data
signal_name (str): the monitor
"""
from lmfit.models import LinearModel, GaussianModel
import matplotlib.pyplot as plt
h = bec.history[hindex]
md = h.metadata["bec"]
# data = h.devices[device_name][signal].read()["value"]
# data_x = h.devices.dcm_pitch.dcm_pitch.read()["value"]
# data_y = h.devices.lu_bpmsum.lu_bpmsum.read()["value"]
data_x = h.devices[device_name][device_name].read()["value"]
data_x = data_x / 1000.0
d0 = 2 * 5.43102 * (1 - 2.4e-4 * 1) / np.sqrt(3.0)
en_vec = 12.39842 / d0 / np.sin(data_x)
data_y = h.devices[signal_name][signal_name].read()["value"]
plt.ion()
plt.figure()
# plt.plot(data_x,data_y,'*')
plt.plot(en_vec, data_y)
plt.title("Energy Scan")
plt.xlabel("E/keV")
plt.ylabel("LU BPM Sum Signal / AU")
plt.show()
dirname = "/sls/x10sa/config/commissioning/Data/"
# writing output to simple data file for later analysis:
combined = np.column_stack((en_vec, data_y))
filename = dirname + "EnScan" + str(hindex) + ".txt"
with open(filename, "w") as f:
np.savetxt(f, combined, delimiter=",", fmt="%5f")
return
#####################
### readmonitors()
#####################
def getdiodepos(diode="i1"):
"""
check if diode is in position
i1 = in KBox, below scinti
i2 = below Eiger
"""
diode_in = 1
dpos = dev.diag_y.user_parameter["i1"] # 44 # mm
measdev = dev.diag_y
diodepos_rb = measdev.user_readback.get()
if abs(dpos - diodepos_rb) > 0.1:
print("Diode not in, please move")
diode_in = 0 # sys.exit(0)
if diode == "i2":
detpos_diode = 214.5
detpos_rb = dev.det_z.user_readback.get()
print("using Det diode requires det up and close and beamstop out")
if abs(detpos_diode - detpos_rb) > 0.5:
print("Diode not in, please move Det up")
diode_in = 0 # sys.exit(0)
return diode_in
def read_mon():
from datetime import datetime
e = getenergy()
fesum = dev.fe_bpmsum.read()["fe_bpmsum"]["value"]
lusum = dev.lu_bpmsum.read()["lu_bpmsum"]["value"]
sssum = dev.ss_bpmsum.read()["ss_bpmsum"]["value"]
# Mono
bragg = dev.dcm_bragg.read()["dcm_bragg"]["value"]
pitch = dev.dcm_pitch.read()["dcm_pitch"]["value"]
perp = dev.dcm_perp.read()["dcm_perp"]["value"]
fpitch = dev.dcm_fpitch.read()["dcm_fpitch"]["value"]
froll = dev.dcm_froll.read()["dcm_froll"]["value"]
gap = dev.id_gap.read()["id_gap"]["value"]
# KB V
vbu = dev.vfm_bu.read()["vfm_bu"]["value"]
vbd = dev.vfm_bd.read()["vfm_bd"]["value"]
vbpitch = dev.vfm_pitch.read()["vfm_pitch"]["value"]
vbyaw = dev.vfm_yaw.read()["vfm_yaw"]["value"]
vbroll = dev.vfm_roll.read()["vfm_roll"]["value"]
vblat = dev.vfm_lat.read()["vfm_lat"]["value"]
vbvert = dev.vfm_vert.read()["vfm_vert"]["value"]
# KB H
hbu = dev.hfm_bu.read()["hfm_bu"]["value"]
hbd = dev.hfm_bd.read()["hfm_bd"]["value"]
hbpitch = dev.hfm_pitch.read()["hfm_pitch"]["value"]
hbyaw = dev.hfm_yaw.read()["hfm_yaw"]["value"]
hbroll = dev.hfm_roll.read()["hfm_roll"]["value"]
hblat = dev.hfm_lat.read()["hfm_lat"]["value"]
hbvert = dev.hfm_vert.read()["hfm_vert"]["value"]
## FE slits centre and size
fe_sx_cen = dev.fe_sxcen.read()["fe_sxcen"]["value"]
fe_sx_size = dev.fe_sxsize.read()["fe_sxsize"]["value"]
fe_sy_cen = dev.fe_sycen.read()["fe_sycen"]["value"]
fe_sy_size = dev.fe_sysize.read()["fe_sysize"]["value"]
## BSF slits centre and size
bsf_xcen = dev.bsf_sl_xcen.read()["bsf_sl_xcen"]["value"]
bsf_xsize = dev.bsf_sl_xsize.read()["bsf_sl_xsize"]["value"]
bsf_ycen = dev.bsf_sl_ycen.read()["bsf_sl_ycen"]["value"]
bsf_ysize = dev.bsf_sl_ysize.read()["bsf_sl_ysize"]["value"]
## SS slits centre and size
ss_sl_xcen = dev.ss_sl_xcen.read()["ss_sl_xcen"]["value"]
ss_sl_xsize = dev.ss_sl_xsize.read()["ss_sl_xsize"]["value"]
ss_sl_ycen = dev.ss_sl_ycen.read()["ss_sl_ycen"]["value"]
ss_sl_ysize = dev.ss_sl_ysize.read()["ss_sl_ysize"]["value"]
## BCU slits centre and size
bcu_sl_xcen = dev.bcu_sl_xcen.read()["bcu_sl_xcen"]["value"]
bcu_sl_xsize = dev.bcu_sl_xsize.read()["bcu_sl_xsize"]["value"]
bcu_sl_ycen = dev.bcu_sl_ycen.read()["bcu_sl_ycen"]["value"]
bcu_sl_ysize = dev.bcu_sl_ysize.read()["bcu_sl_ysize"]["value"]
## move in screen in SS chamber and get size and position
## move out again
# umv(dev.samcam_xmot, 1)
## move in screen at sample position and get size and position
# move up diode and get reading
## move out again
## move in x-ray eye below detector and get size and position
# read out DetDiode and calc flux
## move out again
xeye_x_pos = dev.xeye_x.read()["xeye_x"]["value"]
bcusum = 0
i1signal = 0
if xeye_x_pos < 0:
bcusum = dev.bcu_bpmsum.read()["bcu_bpmsum"]["value"]
i1signal = dev.i1.read()["i1"]["value"]
print("Energy = ", e, " keV")
print(f"fesum,lusum,sssum,bcusum,i1signal = {fesum,lusum,sssum,bcusum,i1signal}")
print(f"bragg, pitch, perp, fpitch, froll, gap = {bragg, pitch,perp, fpitch, froll, gap}")
# return e, fesum,lusum,sssum,bcusum,i1signal
print("KB VERT")
print(
f"vbu, vbd, vbpitch,vbyaw,vbroll,vblat,vbvert = {vbu, vbd, vbpitch,vbyaw,vbroll,vblat,vbvert}"
)
# return vbu, vbd, vbpitch,vbyaw,vbroll,vblat,vbvert
print("KB HORI")
print(
f"hbu, hbd, hbpitch,hbyaw,hbroll,hblat,hbvert = {hbu, hbd, hbpitch,hbyaw,hbroll,hblat,hbvert}"
)
# return hbu, hbd, hbpitch,hbyaw,hbroll,hblat,hbvert
## dump status in CSV
timestamp = datetime.now().strftime("%Y%m%d_%H%M%S")
dirname = "/sls/x10sa/config/commissioning/Data/"
filename = dirname + f"BLstatus_{timestamp}.txt"
with open(filename, "w") as f:
combined = np.column_stack((e, gap))
f.write("En and gap\n")
np.savetxt(f, combined, delimiter=",", fmt="%5f")
f.write("bragg, pitch,perp, fpitch, froll\n")
combined = np.column_stack((bragg, pitch, perp, fpitch, froll))
np.savetxt(f, combined, delimiter=",", fmt="%5f")
f.write("vbu, vbd, vbpitch,vbyaw,vbroll,vblat,vbvert\n")
combined = np.column_stack((vbu, vbd, vbpitch, vbyaw, vbroll, vblat, vbvert))
np.savetxt(f, combined, delimiter=",", fmt="%5f")
f.write("hbu, hbd, hbpitch,hbyaw,hbroll,hblat,hbvert\n")
combined = np.column_stack((hbu, hbd, hbpitch, hbyaw, hbroll, hblat, hbvert))
np.savetxt(f, combined, delimiter=",", fmt="%5f")
f.write("FE slits\n")
combined = np.column_stack((fe_sx_cen, fe_sx_size, fe_sy_cen, fe_sy_size))
np.savetxt(f, combined, delimiter=",", fmt="%5f")
f.write("BSF slits\n")
combined = np.column_stack((bsf_sl_xcen, bsf_sl_xsize, bsf_sl_ycen, bsf_sl_ysize))
np.savetxt(f, combined, delimiter=",", fmt="%5f")
f.write("SS slits\n")
combined = np.column_stack((ss_sl_xcen, ss_sl_xsize, ss_sl_ycen, ss_sl_ysize))
np.savetxt(f, combined, delimiter=",", fmt="%5f")
f.write("BCU slits\n")
combined = np.column_stack((bcu_sl_xcen, bcu_sl_xsize, bcu_sl_ycen, bcu_sl_ysize))
np.savetxt(f, combined, delimiter=",", fmt="%5f")
f.write("fesum,lusum,sssum,bcusum,i1signal\n")
combined = np.column_stack((fesum, lusum, sssum, bcusum, i1signal))
np.savetxt(f, combined, delimiter=",", fmt="%5f")
return
#####################
###long bragg scan
#####################
def longscan():
# scan the entire energy range from 5 keV (bragg angle = 406.512 mrad) to 30 keV (bragg angle = 65.949 mrad)
# "snakewise"
import time
umv(dev.id_gap, 4.5)
time.sleep(0.2)
s = scans.line_scan(dev.dcm_bragg, 406.5, 65.9, steps=1500, exp_time=0.05, relative=False)
time.sleep(2)
umv(dev.id_gap, 5.0)
time.sleep(0.2)
s = scans.line_scan(dev.dcm_bragg, 65.9, 406.5, steps=1500, exp_time=0.05, relative=False)
time.sleep(2)
umv(dev.id_gap, 5.5)
time.sleep(0.2)
s = scans.line_scan(dev.dcm_bragg, 406.5, 65.9, steps=1500, exp_time=0.05, relative=False)
time.sleep(2)
umv(dev.id_gap, 6.0)
time.sleep(0.2)
s = scans.line_scan(dev.dcm_bragg, 65.9, 406.5, steps=1500, exp_time=0.05, relative=False)
time.sleep(2)
umv(dev.id_gap, 6.5)
time.sleep(0.2)
s = scans.line_scan(dev.dcm_bragg, 406.5, 65.9, steps=1500, exp_time=0.05, relative=False)
#####################
### colliscans
#####################
def colliscan(direction: str, range=0.3, nsteps=30, stime=0.5, centre=1):
"""
Colli scan, and possible centre
and plot all the results.
Args:
direction (str): the direction [x/y]
range: default = 0.3
nsteps: default = 30
settling time: default = 0.5
move to centre or not [0/1]
Example: colliscan("h", centre=1)
"""
import sys
dock_area = bec.gui.new()
wr = dock_area.new().new(bec.gui.available_widgets.Waveform)
# check if i1 DIODE is IN
# if not, aks to be moved
diodeinpos = dev.diag_y.user_parameter["i1"] # 44 # mm
colli_up = dev.coll_y.user_parameter["in"] # 41.5 # mm
measdev = dev.diag_y
diodepos_rb = measdev.user_readback.get()
if abs(diodeinpos - diodepos_rb) > 0.1:
print("Diode not in, please move")
sys.exit(0)
signal_name = "i1"
dx = range # scan range
if direction in ("x", "h"):
mot = dev.coll_x
else:
mot = dev.coll_y
# mot_name = "coll_y"
motname = mot.name
# note current motor positions
curr_rbx = dev.coll_x.user_readback.get()
curr_rby = dev.coll_y.user_readback.get()
curr_rb = mot.user_readback.get()
print(f"Motor position before scan: {curr_rb}")
if abs(colli_up - curr_rby) > 1.5:
print(f"Colli vert pos = {curr_rb}, in-pos = {colli_up}")
print("Colli not up, please move")
sys.exit(0)
s = scans.line_scan(mot, -dx, dx, steps=nsteps, settling_time=stime, relative=True)
h = s.scan.data # this is a ScanDataContainer and can get the MetaData
# md = h.metadata["bec"]
data_x = h.devices[motname][motname].read()["value"]
# data_y = h.devices[signal_name][signal_name].read()["value"] # does not work for signal
# data_y = h.devices.signal_name.signal_name.read() # timestamp AND value , dir
data_y = h.devices.i1.i1.read()["value"]
cen, peak = fit_plot(data_x, data_y)
ans = "n"
if abs(cen - curr_rb) > dx:
print("Fit failed, not moving")
ans = input("Going to peak value? y/n ")
if ans == "y":
umv(mot, peak)
else:
if centre:
print(f"Moving to new colli position of {cen}")
umv(mot, cen)
else:
print(f"Scan finished, going back to previous position")
umv(mot, curr_rb)
return
#####################
### slitscans
#####################
def slitscan(device_location: str, direction: str, range: 1, nsteps=50, centre=0):
"""
Slit scan, and possible centre
and plot all the results.
Args:
device_location (str): the slits for which to get the monitoring data
direction (str): the direction [h/v]
nsteps: default = 50
move to centre or not [0/1]
Example: slitscan("fe","h", 30, centre=1)
"""
from lmfit.models import LinearModel, GaussianModel
import matplotlib.pyplot as plt
time = 0.1
dx = range # scan range
# FE slits ================================
if device_location in ["fe"]:
default_h = 1.0 ## 1.4 or 0.4 or ???
default_v = 1.0 ## 1.2
det = dev.lu_bpmsum
if direction == "h":
mot = dev.fe_sxcen
size = dev.fe_sxsize
s_closed = 0.1
s_open = default_h
else:
mot = dev.fe_sycen
size = dev.fe_sysize
s_closed = 0.1
s_open = default_v
# BSF slits ================================
if device_location in ["bsf", "s1"]:
default_h = 3.0 # 8125
default_v = 3.0 # ??? close more ???? # 8149.8
det = dev.lu_bpmsum
# det = dev.ss_bpmsum or #det = dev.bcu_bpmsum would also work
if direction == "h":
mot = dev.bsf_sl_xcen
size = dev.bsf_sl_xsize
s_closed = 0.1
s_open = default_h
else:
mot = dev.bsf_sl_ycen
size = dev.bsf_sl_ysize
s_closed = 0.1
s_open = default_v
# SS slits ================================
if device_location in ["ss", "ss_sl", "ss"]:
default_h = 6.0 # ???
default_v = 5.0 # ??? close more ???? # 8149.8
det = dev.bcu_bpmsum
if direction == "h":
mot = dev.ss_sl_xcen
size = dev.ss_sl_xsize
s_closed = 0.1
s_open = default_h
else:
mot = dev.ss_sl_ycen
size = dev.ss_sl_ysize
s_closed = 0.1
s_open = default_v
# BCU slits ================================
if device_location in ["bcu", "s3", "eb"]:
default_h = 2.0 # ???
default_v = 2.0 # ??? close more ???? # 8149.8
# change to i0 later ??
det = dev.i1
dposm = dev.diag_y.user_parameter["i1"]
# dpos0 = dev.scin_y.user_readback.get()
# if abs(dposm - dpos0) > 1:
# print("moving diode i1 in")
# umv(dev.scin_y, dposm)
if direction == "h":
mot = dev.bcu_sl_xcen
size = dev.bcu_sl_xsize
s_closed = 3.0 ## very large, else does not work !
s_open = default_h
else:
mot = dev.bcu_sl_ycen
size = dev.bcu_sl_ysize
s_closed = 3.0 ## very large !
s_open = default_v
# ===== ================================
print(
f"scanning parameters are: {device_location}, direction= {direction}, range={-range, range}, steps={nsteps}, center it? {centre}"
)
ans = input("All good, start the scan? (y/n): ").strip().lower()
if ans not in ("yes", "y"):
print("Ok, cancelled")
return
dock_area = bec.gui.new()
wr = dock_area.new().new(bec.gui.available_widgets.Waveform)
pos0 = mot.user_readback.get()
siz0 = size.user_readback.get()
umv(size, s_closed)
s = scans.line_scan(mot, -dx, dx, steps=nsteps, exp_time=time, relative=True)
wr.plot(x_name=mot.name, y_name=det.name)
wr.x_label = mot.name
wr.y_label = det.name
# this gives funny errors: ARRAYS differently read?
# data_x = s.scan.data.devices[mot.name][mot.name].read()["value"]
# data_x = s.scan.live_data["fe_sxcen"]["fe_sxcen"].val
# data_y = s.scan.data.devices[det.name][det.name].read()["value"]
# FE slits ================================
data_y = s.scan.live_data.lu_bpmsum.lu_bpmsum.val
if mot.name == "fe_sxcen":
data_x = s.scan.live_data.fe_sxcen.fe_sxcen.val
if mot.name == "fe_sycen":
data_x = s.scan.live_data.fe_sycen.fe_sycen.val
# BSF slits ================================
if mot.name == "bsf_sl_xcen":
data_x = s.scan.live_data.bsf_sl_xcen.bsf_sl_xcen.val
if mot.name == "bsf_ycen":
data_x = s.scan.live_data.bsf_sl_ycen.bsf_sl_ycen.val
# SS slits ================================
if mot.name == "ss_sl_xcen":
data_x = s.scan.live_data.ss_sl_xcen.ss_sl_xcen.val
data_y = s.scan.live_data.bcu_bpmsum.bcu_bpmsum.val
if mot.name == "ss_sl_ycen":
data_x = s.scan.live_data.ss_sl_ycen.ss_sl_ycen.val
data_y = s.scan.live_data.bcu_bpmsum.bcu_bpmsum.val
# BCU slits ================================
if mot.name == "bcu_sl_xcen":
data_x = s.scan.live_data.bcu_sl_xcen.bcu_sl_xcen.val
data_y = s.scan.live_data.i1.i1.val
if mot.name == "bcu_sl_ycen":
data_x = s.scan.live_data.bcu_sl_ycen.bcu_sl_ycen.val
data_y = s.scan.live_data.i1.i1.val
# change to i0 later ??
wr.plot(data_x, data_y)
# fit Gaussian + lin BG
peak = GaussianModel()
background = LinearModel()
model = peak + background
maxy = max(data_y)
indmax = np.argmax(data_y)
xm = data_x[indmax]
p = model.make_params(amplitude=maxy, center=xm)
g = model.fit(data_y, p, x=data_x)
print(f'Center: {g.params["center"].value:.5f} FWHM: {g.params["fwhm"].value:.5f}')
plt.ion()
plt.figure(2)
plt.plot(data_x, data_y, "*")
plt.plot(data_x, g.best_fit)
plt.show()
res_center = g.params["center"].value
if centre and abs(res_center) < 1.5:
print(f"moving slit centre from {pos0} to {res_center}")
umv(mot, res_center)
else:
print("not changing the slit position")
print(f"opening slit to {s_open}")
umv(size, s_open)
return
#####################
### KB focusing
#####################
## note acceptance of KB;
## length = 400 mm
## about 3 mrad pitch --> roughly 1.5 mm in both dir
##
def kbfocus(sizex, sizey):
#
# note: bending is in neg direction, and the KB mirrors have a hysteresis
# better to bend from ONE direction only, always unbend/bend for more reproducible results
#
import math
import time
en = getenergy()
print(f"Energy is {en} keV")
print("Currently only 2 sizes supported, small approx.(2 x 2.7) and medium approx.(40 x 40)")
## is the pitch ok ?
vpitch = 2.695
hpitch = 2.997
accepted_tol = 5e-3
vpitch_rb = dev.vfm_pitch.user_readback.get()
hpitch_rb = dev.hfm_pitch.user_readback.get()
if not math.isclose(vpitch_rb, vpitch, rel_tol=accepted_tol):
print("warning! vertical pitch may be off")
if not math.isclose(hpitch_rb, hpitch, rel_tol=accepted_tol):
print("warning! horizontal pitch may be off")
##### one good size:
#### careful, the scinti position !
vfbend_u = -0.8212
vfbend_d = -1.3250
hfbend_u = -2.7270
hfbend_d = -5.3480
# if sizey < 5:
# vfbend_u = -0.9895 # -1.029
# vfbend_d = -1.4050 # -1.545
# else:
# vfbend_u = -0.9793
# vfbend_d = -1.4950
#
# if sizex < 5:
# hfbend_u = -2.887
# hfbend_d = -5.550
# else:
# hfbend_u = -2.7378
# hfbend_d = -5.4000
umv(dev.vfm_bu, vfbend_u + 0.01)
umv(dev.vfm_bd, vfbend_d + 0.01)
time.sleep(0.1)
umv(dev.vfm_bu, vfbend_u)
umv(dev.vfm_bd, vfbend_d)
umv(dev.hfm_bu, hfbend_u + 0.01)
umv(dev.hfm_bd, hfbend_d + 0.01)
time.sleep(0.1)
umv(dev.hfm_bu, hfbend_u)
umv(dev.hfm_bd, hfbend_d)
return
##### ======================
##### goto several energies with about the right gap opening
def goto7():
perp_dist = 0.0389
umv(dev.dcm_perp, perp_dist)
gap = 6.9 # 3rd harm
umv(dev.id_gap, gap)
sete(7)
def goto20():
perp_dist = 0.0389
umv(dev.dcm_perp, perp_dist)
gap = 4.736 # 13th harm
umv(dev.id_gap, gap)
sete(20)
def goto27():
perp_dist = 0.0389
umv(dev.dcm_perp, perp_dist)
gap = 4.86 # 17th harm
umv(dev.id_gap, gap)
sete(27)
def goto30():
perp_dist = 0.0389
umv(dev.dcm_perp, perp_dist)
gap = 4.842 # 19th harm
umv(dev.id_gap, gap)
sete(30)
#########################
### Status of the BEC
#########################
def bstatus():
#
dock_area = bec.gui.new()
dbrowser = dock_area.new("device_browser").new(bec.gui.available_widgets.DeviceBrowser)
dock_area.new("queue").new(bec.gui.available_widgets.BECQueue)
# queue = dock_area.queue.BECQueue # give it a name
# text_box = dock_area.new().new(widget=bec.gui.available_widgets.TextBox)
# text_box.set_plain_text("Hello, World!")
# text_box.set_html_text("<h1>Welcome to BEC Widgets</h1><p>This is an example of displaying <strong>HTML</strong> text.</p>")
###########################
### move in the X-ray Eye
###########################
def detxeye_in():
"""
moves in the Xray eye and the Diode below detector
"""
#
setfoc = 1.0
setzoom = 0.1
setxrpos = 2.0
setdetzpos = 260 # 200 better?
setdetypos = 214.5
detzpos = dev.det_z.user_readback.get()
detypos = dev.det_y.user_readback.get()
print(f"detector position now: z = {detzpos}, y = {detypos}")
ans = "n"
ans = input("Is it ok to move detector in and up? y/n ")
if ans == "y":
print("moving the detector up (only while detz GT 300 mm) and in")
if detzpos > 300:
umv(dev.det_y, setdetypos)
umv(dev.det_z, setdetzpos)
else:
print("ok, doing nothing")
return
print("moving X-ray eye below Eiger in")
xrpos = dev.det_xi_x.user_readback.get()
if abs(setxrpos - xrpos) > 1:
umv(dev.det_xi_x, setxrpos)
zoompos = dev.det_xi_zoom.user_readback.get()
if abs(setzoom - zoompos) > 1:
umv(dev.det_xi_zoom, setzoom)
focpos = dev.det_xi_focus.user_readback.get()
if abs(setfoc - focpos) > 1:
umv(dev.det_xi_focus, setfoc)
def detxeye_out():
"""
moves det back and down
"""
setdetzpos = 500
setdetypos = 0
umv(dev.det_z, setdetzpos)
detzpos = dev.det_z.user_readback.get()
if detzpos > 300:
umv(dev.det_y, setdetypos)
###########################
### beamsize on scinti
###########################
def measure_samcam(zoom=1000):
scinti_inpos = 38.6 # mm
sc_rb = dev.diag_y.user_readback.get()
if abs(scinti_inpos - sc_rb) > 0.3:
print("Scinti not in, please move")
sys.exit(0)
auto_exposure(cam="samcam", target=200)
a = dev.samcam_xsig.read()["samcam_xsig"]["value"]
b = dev.samcam_ysig.read()["samcam_ysig"]["value"]
sx = calc_scam_microns(a, zoom=zoom) * 2.35
sy = calc_scam_microns(b, zoom=zoom) * 2.35
print(f"FWHM in um : {sx, sy}")
return sx, sy
def measure_sscam():
x_inpos = dev.ss_xi_x.user_parameter["in"] # 7 # mm
px2mum = 20
scpos_rb = dev.ss_xi_x.user_readback.get()
if abs(x_inpos - scpos_rb) > 0.3:
print("Scinti not in, please move")
sys.exit(0)
auto_exposure(cam="sscam", target=200)
a = dev.ss_xicam_xsig.read()["ss_xicam_xsig"]["value"]
b = dev.ss_xicam_ysig.read()["ss_xicam_ysig"]["value"]
sx = a * px2mum * 2.35
sy = b * px2mum * 2.35
print(f"FWHM at SS cam in um : {sx, sy}")
return sx, sy
######### special purpose programs ###############
##################################################
# EXAMPLE: knife_scan.py
def knife_edge(dir="hor", range=0.05, steps=100):
mot = dev.gon_x
if dir == "vert":
mot = dev.gon_y
s = scans.line_scan(dev.gon_x, -range, range, steps=steps, exp_time=1, relative=True)
return
##############
def smoothplotgrad(dx, dy, reverse=0):
"""
for example from the knife edge scans
if reverse == 1 then mirrored for plot and fit
USE as: smoothplotgrad(dx,dy, reverse = 0)
"""
from scipy.ndimage import gaussian_filter1d
import matplotlib.pyplot as plt
gf = gaussian_filter1d(dy, 1)
plt.ion()
plt.figure()
plt.plot(dx, dy, "*")
plt.plot(dx, gf)
plt.figure()
grf = np.gradient(gf)
if reverse:
grf = -grf
plt.plot(dx, grf)
g, xm = justfit(dx, grf, model="gauss", ibg=0)
plt.plot(dx, g.best_fit)
return
##################################################################################
### enery scan program around a harmonic
def scan_eg(erange, nsteps=50, fit=True):
"""
Scan an energy range around a certain harmonic
scan_eg(erange, steps= 50, fit= True) # default
scan_eg(erange, steps= 50, fit= False) # no fit wanted
eg, erange = 12.4 *0.02 for a width of 2 % around 12.4 keV
"""
import time
from datetime import datetime
from lmfit.models import LinearModel, GaussianModel, VoigtModel
import matplotlib.pyplot as plt
## defaults
exptime = 1
e = getenergy()
print(f"Current energy is: {e} keV")
e_start = e - erange
e_end = e + erange
en_resol = (e_start - e_end) / nsteps
print(f"energy resolution is: {en_resol} keV")
mot_scan = dev.dcm_bragg
mot_pitch = dev.dcm_pitch
mot_perp = dev.dcm_perp
print(
f"scanning parameters are: motor = {mot_scan.name}, range={e-erange, e+erange}, steps={nsteps}"
)
ans = input("All good, start the scan? (y/n): ").strip().lower()
if ans not in ("yes", "y"):
print("Ok, cancelled")
return
mot_scan_rb = mot_scan.user_readback.get()
mot_pitch_rb = mot_pitch.user_readback.get()
mot_perp_rb = mot_perp.user_readback.get()
det = dev.lu_bpmsum
print(f"Position of bragg angle is {mot_scan_rb} mrad")
print(f"Position of pitch angle is {mot_pitch_rb} mrad")
print(f"Perp distance is {mot_perp_rb} mm")
a_start = angle(e_start, "ln") * 1000
a_end = angle(e_end, "ln") * 1000.0
print(f"Scanning Bragg from {a_start} to {a_end} mrad")
s = scans.line_scan(mot_scan, a_start, a_end, steps=nsteps, exp_time=exptime, relative=False)
## plot and fit the scan
bragg_data = (
s.scan.live_data.dcm_bragg.dcm_bragg.val
) # is a list and needs to be converted to array
bpm_data = s.scan.live_data.lu_bpmsum.lu_bpmsum.val
if fit:
print("Bragg Angle")
bcen, xm = fit_plot(bragg_data, bpm_data, model="voigt", fitrange=1)
print("")
br_rad = np.array(bragg_data) / 1000.0
d0 = 2 * 5.43102 * (1 - 2.4e-4 * 1) / np.sqrt(3.0)
energy_data = 12.39842 / d0 / np.sin(br_rad)
print("Energy")
ecen, xm = fit_plot(energy_data, bpm_data, model="voigt", fitclick=1)
## save data to file ?
ans = "n"
ans = input("Store data in csv file? y/n ")
if ans == "y":
timestamp = datetime.now().strftime("%Y%m%d_%H%M%S")
data_x = bragg_data
data_e = energy_data
data_y = bpm_data
dirname = "/home/gac-x10sa/Data/"
# writing output to simple data file for later analysis:
combined = np.column_stack((data_x, data_e, data_y))
tname = f"Scan_{timestamp}.txt"
filename = dirname + tname
with open(filename, "w") as f:
np.savetxt(f, combined, delimiter=",", fmt="%5f")
sete(e) # go back
return
###########################
### compute a signal
###########################
# see in config file
################################################
### open window/doch for long gap scan
################################################
#
### Note: all deprecated !!!!
def scan_window(wname="Scan", fit=True):
dock_area = bec.gui.new()
# Add a new dock with a Waveform to the BECDockArea
nam = "waveform_dock_"+wname
dock_area.new(name=nam, widget="Waveform")
# dock_area.panels
# dock_area.panel_list
plt1 = dock_area.panels[nam]
# Add signals to the WaveformWidget
plt1.plot(device_x='id_gap', device_y='bpm')
dock2 = dock_area.new(name="motor_dock", widget="MotorMap",relative_to=nam, position="right")
###do stuff
### if done, remove
dock2.remove()