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)