Files
sf-op/script/RFscan/phase_scan_caqtdm.py
2026-02-09 15:06:52 +01:00

195 lines
9.1 KiB
Python

import ch.psi.pshell.epics.Positioner as Positioner
import ch.psi.pshell.epics.ChannelDouble as ChannelDouble
import_py("CPython/hfitoff", "hfitoff")
dry_run = False
do_elog = True
#station = "SINSB02" # define param locally
#bpm_ch = "SINBC02-DBPM140"
#bpm_field = "X1"
run("RFscan/phase_scan_data.py")
if get_exec_pars().args: # args is set by callin process (Qt panel)
station = args[0]
bpm_ch = args[1]
bpm_field = args[2]
if dry_run:
station = "STEST01"
bpm_ch = "SINBC02-DBPM140"
bpm_field = "X1-SIMU"
bpm_val = ChannelDouble("BPM position", bpm_ch + ":" + bpm_field)
amplt = ChannelDouble(station + " Amplitude", station + "-RSYS:GET-VSUM-AMPLT")
power = ChannelDouble(station + " Klystron Power", station + "-RSYS:GET-KLY-POWER-GATED")
bpm_val.initialize()
amplt.initialize()
power.initialize()
phase = Positioner("Phase", station + "-RSYS:SET-VSUM-PHASE", station + "-RSYS:GET-VSUM-PHASE")
phase.config.minValue = -90.0 # cannot be zero, we need margin for scanning around 0 deg
phase.config.maxValue = 360.0
phase.config.precision = 4 # digits beyond this are ignored
phase.config.deadband = 0.3 # dead-band for comparing set and get values
phase.config.rotation = True # forced to 0..360 deg for comparison
phase.initialize()
phase0 = phase.read()
start = caget(station + "-RSYS:SET-SCAN-START")
end = caget(station + "-RSYS:SET-SCAN-STOP")
step = caget(station + "-RSYS:SET-SCAN-STEP")
lat = caget(station + "-RSYS:SET-SCAN-WAIT-TIME")
nb = caget(station + "-RSYS:SET-NUM-AVERAGE")
disp = caget(bpm_ch + ":DISPERSION-OP")
p0 = caget(bpm_list[bpm_ch]["mbnd"] + ":P-READ")
bpm_averager = create_averager(bpm_val, nb, lat)
energy0 = p0 - 0.511
A, B = energy0 / (disp * 1000), energy0 # param to convert bpm mm to MeV
# pause scan when beam is down
def is_beam_ok():
mps_ma = int(caget("SIN-EMED-DTSA001:LEVEL1_MA.RVAL"))
mps_ar = int(caget("SIN-EMED-DTSA002:LEVEL1_AR.RVAL"))
mps_at = int(caget("SIN-EMED-DTSA003:LEVEL1_AT.RVAL"))
beam_ok = (mps_ma + mps_ar + mps_at == 3)
return beam_ok
def wait_beam():
if not is_beam_ok():
print("Waiting for beam...")
while not is_beam_ok():
time.sleep(0.5)
print("Beam ok")
def before(rec):
wait_beam()
# check for beam presence and update the plot dynamically
arr_phase, arr_energy = [],[]
def after(rec):
if not is_beam_ok():
print("Beam is down, invalidating record...")
rec.invalidate()
global A, B
arr_phase.append(rec.positions[0])
arr_energy.append(A * rec.readables[0].mean + B)
caput(station + "-RSYS:GET-PHASE-ARRAY", to_array(arr_phase, 'd'))
caput(station + "-RSYS:GET-ENERGY-ARRAY", to_array(arr_energy,'d'))
caput("SF-PHASE-GLOBAL:GET-PHASE-ARRAY", to_array(arr_phase, 'd'))
caput("SF-PHASE-GLOBAL:GET-ENERGY-ARRAY", to_array(arr_energy,'d'))
# scan and plot
try:
feedback_state_init = {}
for feedback_channel in station_list[station]["feedback"]:
feedback_state = caget(feedback_channel)
feedback_state_init[feedback_channel] = feedback_state
caput(feedback_channel, 0)
phase_tol_init = caget(station + "-RSYS:SET-PHASE-TOL")
if phase_tol_init < 10.0 :
caput(station + "-RSYS:SET-PHASE-TOL", 10.0)
nominal_beam_phase = 270.0 if station == "SINXB01" else 90.0
caput(station + "-RSYS:SET-BEAM-PHASE", nominal_beam_phase)
time.sleep(1.0)
ph_crest0 = phase.read()
phase.write(start)
time.sleep(2.0)
caput(station + "-RSYS:PHASE-SCAN-MESSAGE", "scanning " + station)
caput(station + "-RSYS:GET-ONCREST-VSUM-PHASE", float('nan'))
caput(station + "-RSYS:GET-ONCREST-VSUM-AMPLT", float('nan'))
caput(station + "-RSYS:GET-ONCREST-E-GAIN", float('nan'))
caput(station + "-RSYS:GET-ONCREST-KLY-POWER", float('nan'))
caput(station + "-RSYS:GET-FIT-PHASE-ARRAY", to_array([0.0],'d'))
caput(station + "-RSYS:GET-FIT-ENERGY-ARRAY", to_array([0.0],'d'))
caput("SF-PHASE-GLOBAL:GET-FIT-PHASE-ARRAY", to_array([0.0],'d'))
caput("SF-PHASE-GLOBAL:GET-FIT-ENERGY-ARRAY", to_array([0.0],'d'))
r = lscan(phase, bpm_averager, start, end, step, latency=lat, before_read = before, after_read = after)
energy = [A * val.mean + B for val in r.getReadable(0)] # convert bpm mm to MeV
rf_phase = r.getPositions(0)
if start < 0:
rf_phase = [((ph + 90) % 360) -90 for ph in rf_phase ] # force phase values in -90..270 range
caput(station + "-RSYS:GET-ENERGY-ARRAY", to_array(energy, 'd'))
caput(station + "-RSYS:GET-PHASE-ARRAY", to_array(rf_phase,'d'))
caput("SF-PHASE-GLOBAL:GET-ENERGY-ARRAY", to_array(energy, 'd'))
caput("SF-PHASE-GLOBAL:GET-PHASE-ARRAY", to_array(rf_phase,'d'))
try:
fit_amplitude, fit_phase_deg, fit_offset, ph_crest, fit_x, fit_y = hfitoff(energy, rf_phase)
except:
raise Exception("Fit failure")
plt = plot(None, name="phase scan")[0]
plt.getSeries(0).setData(to_array(rf_phase,'d'), to_array(energy,'d'))
plt.getSeries(0).setPointSize(6)
plt.getSeries(0).setLinesVisible(False)
plt.addSeries(LinePlotSeries("fit"))
plt.getSeries(1).setData(fit_x, fit_y)
plt.getSeries(1).setPointsVisible(False)
plt.setLegendVisible(True)
ph_crest = (ph_crest + 180) % 360 if station == "SINXB01" else ph_crest # take "on-crest" as 270 deg for SINXB01
phase.write(ph_crest)
time.sleep(1.0)
nominal_ampl = amplt.read()
nominal_power = power.read()
caput(station + "-RSYS:GET-ONCREST-VSUM-PHASE", ph_crest)
caput(station + "-RSYS:GET-ONCREST-E-GAIN", fit_amplitude)
caput(station + "-RSYS:GET-ONCREST-VSUM-AMPLT", nominal_ampl)
caput(station + "-RSYS:GET-ONCREST-KLY-POWER", nominal_power)
caput(station + "-RSYS:GET-FIT-PHASE-ARRAY", fit_x)
caput(station + "-RSYS:GET-FIT-ENERGY-ARRAY", fit_y)
caput("SF-PHASE-GLOBAL:GET-FIT-PHASE-ARRAY", fit_x)
caput("SF-PHASE-GLOBAL:GET-FIT-ENERGY-ARRAY", fit_y)
now = str(caget("SF-TIME:FULL-DATE"))
now = now[:16] if now[0].isdigit() else ""
message = now + " " + str('{:.2f}'.format(ph_crest)) + " deg"
message = message + " (" + str('{:.2f}'.format(ph_crest - ph_crest0)) + ")"
caput(station + "-RSYS:PHASE-SCAN-MESSAGE", message)
phase_corr = caget(station + "-RSYS:GET-VSUM-PHASE-OFFSET-CORR")
phase_offset = nominal_beam_phase - ph_crest - phase_corr
amplitude_scale = fit_amplitude / nominal_ampl if nominal_ampl != 0 else 0.0
power_scale = nominal_power / fit_amplitude**2 if fit_amplitude != 0 else 0.0
caput(station + "-RSYS:SET-VSUM-PHASE-OFFSET-BASE-CALC", phase_offset)
caput(station + "-RSYS:SET-VSUM-AMPLT-SCALE-CALC", amplitude_scale)
caput(station + "-RSYS:SET-VOLT-POWER-SCALE-CALC", power_scale)
finally:
for feedback_channel in station_list[station]["feedback"]:
feedback_state = feedback_state_init[feedback_channel]
caput(feedback_channel, feedback_state)
caput(station + "-RSYS:SET-PHASE-TOL", phase_tol_init)
phase.write(phase0)
Ph_Beam = caget(station + "-RSYS:SET-BEAM-PHASE")
while Ph_Beam > 360:
Ph_Vsum = Ph_Beam - 360
while Ph_Beam < 0:
Ph_Beam = Ph_Beam + 360
caput(station + "-RSYS:SET-BEAM-PHASE", Ph_Beam)
phase.close()
amplt.close()
power.close()
bpm_val.close()
# Saving metadata
save_dataset ("experiment/Station" , station )
save_dataset ("scan 1/processed/Energy" , energy )
save_dataset ("scan 1/processed/Energy gain" , fit_amplitude )
save_dataset ("scan 1/processed/On-crest VS-phase" , ph_crest )
save_dataset ("scan 1/processed/VS-phase offset" , phase_offset )
save_dataset ("scan 1/processed/Amplitude scale" , amplitude_scale )
save_dataset ("scan 1/processed/Power scale" , power_scale )
save_dataset ("scan 1/processed/Dispersion" , disp )
set_attribute("scan 1/processed/Energy" , "Unit", "MeV" )
set_attribute("scan 1/processed/Energy gain" , "Unit", "MeV" )
set_attribute("scan 1/processed/On-crest VS-phase" , "Unit", "deg" )
set_attribute("scan 1/processed/VS-phase offset" , "Unit", "deg" )
set_attribute("scan 1/processed/Amplitude scale" , "Unit", "MV" )
set_attribute("scan 1/processed/Power scale" , "Unit", "MW/MV^2" )
set_attribute("scan 1/processed/Dispersion" , "Unit", "mm" )
# Elog entry
if do_elog:
title = "Phase scan " + station
log_msg = "Data file: " + get_exec_pars().path + "\n"
log_msg = log_msg + "On-crest VS phase: %0.2f" % ph_crest + " deg \n"
log_msg = log_msg + "Energy gain: %0.3f" % fit_amplitude + " MeV \n"
log_msg = log_msg + "VS-phase offset: %0.2f" % phase_offset + " deg \n"
log_msg = log_msg + "Amplitude scale: %0.3f" % amplitude_scale + " MV \n"
log_msg = log_msg + "Power scale: %0.7f" % power_scale + " MW/MV^2 \n"
log_msg = log_msg + "Dispersion: %0.4f" % disp + " mm"
attachments = get_plot_snapshots(size=(600,400))
elog(title, log_msg, attachments)