From 4899bf0bd87df02756d31729c21283bf454f8fd3 Mon Sep 17 00:00:00 2001 From: gac bernina Date: Sat, 11 Jul 2026 14:08:43 +0200 Subject: [PATCH] timetool calib work and scan correctness --- eco/acquisition/daq_client.py | 84 +++++++++---- eco/acquisition/scan_data.py | 221 +++++++++++++++++++++++++++++++++- eco/timing/timing_diag.py | 123 ++++++++++++++++--- 3 files changed, 387 insertions(+), 41 deletions(-) diff --git a/eco/acquisition/daq_client.py b/eco/acquisition/daq_client.py index f0a32e1..885b15a 100644 --- a/eco/acquisition/daq_client.py +++ b/eco/acquisition/daq_client.py @@ -1,7 +1,7 @@ import json import pickle import shutil -from threading import Thread +from threading import Thread, Lock, Event import time import traceback import colorama @@ -73,11 +73,23 @@ class Daq(Assembly): self._pgroup = pgroup if type(pulse_id_adj) is str: self.pulse_id = DetectorPvDataStream(pulse_id_adj, name="pulse_id") - self._pid_wo_automonitor = PV( - "SGE-CPCW-85-EVR0:RX-PULSEID", - connection_timeout=0.05, - auto_monitor=False, - ) + # Dedicated, permanently-monitored PV used only to wait for a fresh + # pulse_id in start() without issuing competing CA get requests. + # Kept separate from self.pulse_id._pv (auto_monitor=False) so that + # get_current_value() elsewhere is unaffected; both share the same + # underlying CA channel, so this costs nothing extra on the wire. + self._pulse_id_latest = {"value": None, "timestamp": None} + self._pulse_id_latest_lock = Lock() + self._pulse_id_updated = Event() + + def _on_pulse_id_update(value=None, timestamp=None, **kwargs): + with self._pulse_id_latest_lock: + self._pulse_id_latest["value"] = value + self._pulse_id_latest["timestamp"] = timestamp + self._pulse_id_updated.set() + + self._pulse_id_monitor_pv = PV(pulse_id_adj, auto_monitor=True) + self._pulse_id_monitor_pv.add_callback(_on_pulse_id_update) else: self.pulse_id = pulse_id_adj self.running = [] @@ -246,26 +258,50 @@ class Daq(Assembly): ) starttime_local = time.time() - pv = self.pulse_id._pv start_id = None - tvars = None - poll_interval = 0.02 - max_poll_interval = 0.25 - per_call_timeout = 1.0 # headroom for a saturated CA processing thread - while True: - tvars = pv.get_timevars(timeout=per_call_timeout) - if tvars is not None and tvars["timestamp"] >= starttime_local: - start_id = pv.get(use_monitor=False, timeout=per_call_timeout) - if start_id is not None: + if hasattr(self, "_pulse_id_updated"): + # Wait on the dedicated pulse_id monitor's cache instead of issuing + # explicit CA get requests, to avoid adding CA traffic that competes + # with itself/other concurrent CA activity right at scan start. + while True: + with self._pulse_id_latest_lock: + ts = self._pulse_id_latest["timestamp"] + val = self._pulse_id_latest["value"] + if ts is not None and val is not None and ts >= starttime_local: + start_id = int(val) break - if time.time() - starttime_local > self.timeout: - raise TimeoutError( - f"Timeout {self.timeout} s hit while waiting for a valid, up-to-date " - f"pulse_id. timevars: {tvars}; start_id: {start_id}; " - f"starttime of scan step: {starttime_local}" - ) - time.sleep(poll_interval) - poll_interval = min(poll_interval * 1.5, max_poll_interval) + remaining = self.timeout - (time.time() - starttime_local) + if remaining <= 0: + raise TimeoutError( + f"Timeout {self.timeout} s hit while waiting for a valid, " + f"up-to-date pulse_id. last timestamp: {ts}; " + f"starttime of scan step: {starttime_local}" + ) + self._pulse_id_updated.wait(timeout=min(remaining, 0.25)) + self._pulse_id_updated.clear() + else: + # Fallback for a pulse_id_adj configured as a pre-built object rather + # than a PV name string: no dedicated monitor is set up in that case, + # so fall back to an explicit polling wait. + pv = self.pulse_id._pv + tvars = None + poll_interval = 0.02 + max_poll_interval = 0.25 + per_call_timeout = 1.0 # headroom for a saturated CA processing thread + while True: + tvars = pv.get_timevars(timeout=per_call_timeout) + if tvars is not None and tvars["timestamp"] >= starttime_local: + start_id = pv.get(use_monitor=False, timeout=per_call_timeout) + if start_id is not None: + break + if time.time() - starttime_local > self.timeout: + raise TimeoutError( + f"Timeout {self.timeout} s hit while waiting for a valid, up-to-date " + f"pulse_id. timevars: {tvars}; start_id: {start_id}; " + f"starttime of scan step: {starttime_local}" + ) + time.sleep(poll_interval) + poll_interval = min(poll_interval * 1.5, max_poll_interval) acq_pars = { "label": label, diff --git a/eco/acquisition/scan_data.py b/eco/acquisition/scan_data.py index b8d91fd..fd349fb 100644 --- a/eco/acquisition/scan_data.py +++ b/eco/acquisition/scan_data.py @@ -7,6 +7,7 @@ from escape.swissfel import load_dataset_from_scan import json from pandas import DataFrame +from tabulate import tabulate class RunData: @@ -181,11 +182,224 @@ class StatusData: return r +def _walk_status_items(base_obj, node, depth, include_hidden, _seen): + """Recursively collect (name relative to base_obj, item) pairs from + node's status_collection, unfolding nested assemblies while depth + allows (depth<=1 stops unfolding further). Mirrors the traversal used + by Assembly.get_tree/get_display_str, but callable at arbitrary depth + and independent of the "display" selection when include_hidden=True. + """ + rows = [] + display_names = set(node.status_collection.selections.get("display", {}).keys()) + for wref in node.status_collection._list: + item = wref() + if item is None or item is node: + continue + if id(item) in _seen: + continue + try: + local_name = item.alias.get_full_name(base=node) + except Exception: + continue + if not include_hidden and local_name not in display_names: + continue + _seen.add(id(item)) + try: + name = item.alias.get_full_name(base=base_obj) + except Exception: + name = local_name + is_nested = hasattr(item, "status_collection") + if is_nested and depth > 1: + sub_rows = _walk_status_items(base_obj, item, depth - 1, include_hidden, _seen) + if sub_rows: + rows.extend(sub_rows) + continue + rows.append((name, item)) + return rows + + +class RunStatusView: + """Historical, run-scoped counterpart to an Assembly's live display + repr (`get_display_str`/`__repr__`). Returned by RunStatusAccessor, + e.g. `obj.run_status(1234)` or `obj.run_status[1234]` - not meant to + be constructed directly. + """ + + def __init__( + self, obj, run_number, dic, depth=None, include_hidden=False, compare_live=False + ): + self.obj = obj + self.run_number = run_number + self.dic = dic + self.depth = depth + self.include_hidden = include_hidden + self.compare_live = compare_live + + def _rows(self): + obj = self.obj + if not hasattr(obj, "status_collection"): + return [(name, None) for name in self.dic] + if not self.include_hidden and self.depth is None: + # default: exactly the items the live repr would show + items = obj.status_collection.get_list(selection="display") + return [(item.alias.get_full_name(base=obj), item) for item in items] + return _walk_status_items(obj, obj, self.depth or 1, self.include_hidden, set()) + + def to_dict(self): + """Raw stored values for this run, as a flat dict (no display formatting).""" + return dict(self.dic) + + def _build_table(self, tablefmt="simple"): + tab = [] + for name, item in self._rows(): + value = self.dic.get(name) + if value is None: + prefix = name + "." + if any(k.startswith(prefix) for k in self.dic): + value = "\x1b[3mnested - unfold with depth=/include_hidden=\x1b[0m" + else: + value = "" + typechar = "" + unit = "" + description = "" + if item is not None: + if hasattr(item, "status_collection"): + typechar = "↳" + try: + unit = item.unit.get_current_value() + except Exception: + pass + try: + description = item.description.get_current_value() + except Exception: + pass + if self.compare_live: + try: + value = f"{value} (crr: {item.get_current_value()})" + except Exception: + pass + tab.append([name, value, unit, typechar, description]) + if not tab: + return "" + return tabulate(tab, tablefmt=tablefmt, maxcolwidths=[None, 50, None, None, None]) + + def __getitem__(self, key): + return self.dic[key] + + def __repr__(self): + name = self.obj.alias.get_full_name() + header = f"{name} run {self.run_number} status" + if self.compare_live: + header += " [value shown as: run value (crr: )]" + return header + "\n" + self._build_table() + + def _repr_html_(self): + return self._build_table(tablefmt="html") + + +class RunStatusAccessor: + """`obj.run_status` - pick a run to see this object's historical status + as a display view (like its live repr), with options to unfold nested + (`depth=`) and normally-hidden (`include_hidden=`) components, and to + compare against the live value (`compare_live=`). Tab-completable in + IPython/Jupyter via `obj.run_status[]`. Use `.dict(...)` for the + original, raw (optionally multi-run) dictionary/DataFrame output. + """ + + def __init__(self, obj): + self._obj = obj + + def _pgroup_key(self, pgroup): + if pgroup == "auto": + return list(STATUS_DATA.keys())[0] + return pgroup + + def available(self, pgroup="auto"): + return STATUS_DATA[self._pgroup_key(pgroup)].get_available_run_numbers() + + def keys(self, pgroup="auto"): + """Enables dict-key tab completion, e.g. obj.run_status[].""" + return self.available(pgroup=pgroup) + + def dict( + self, + run_number=None, + par_type="status", + force_reload=False, + pgroup="auto", + status_type="status_run_start", + as_dataframe=False, + ): + """Raw values for one or several runs, as {run_number: {name: value}} + (or a DataFrame if as_dataframe=True). This is the original + run_status() return shape.""" + return self._obj._fetch_run_status_dict( + run_number=run_number, + par_type=par_type, + force_reload=force_reload, + pgroup=pgroup, + status_type=status_type, + as_dataframe=as_dataframe, + ) + + def __call__( + self, + run_number=None, + depth=None, + include_hidden=False, + compare_live=False, + par_type="status", + force_reload=False, + pgroup="auto", + status_type="status_run_start", + ): + if run_number is None: + return self.available(pgroup=pgroup) + single = isinstance(run_number, Number) + run_numbers = [run_number] if single else list(run_number) + stat = self._obj._fetch_run_status_dict( + run_number=run_numbers, + par_type=par_type, + force_reload=force_reload, + pgroup=pgroup, + status_type=status_type, + ) + views = { + runno: RunStatusView( + self._obj, + runno, + dic, + depth=depth, + include_hidden=include_hidden, + compare_live=compare_live, + ) + for runno, dic in stat.items() + } + if single: + return list(views.values())[0] + return views + + def __getitem__(self, run_number): + return self(run_number) + + def __repr__(self): + try: + runs = self.available() + avail = f"{len(runs)} runs available ({min(runs)}-{max(runs)})" if runs else "no runs available" + except Exception: + avail = "run list unavailable" + return ( + f"\n" + "Call e.g. .run_status(1234) or .run_status[1234] for a display view " + "(depth=, include_hidden=, compare_live=); .run_status.dict(...) for raw values." + ) + + def run_status_convenience(Obj): # if not hasattr(Obj, "alias"): # return Obj - def run_status( + def _fetch_run_status_dict( self, run_number=None, par_type="status", @@ -223,7 +437,8 @@ def run_status_convenience(Obj): return stat - Obj.run_status = run_status + Obj._fetch_run_status_dict = _fetch_run_status_dict + Obj.run_status = property(lambda self: RunStatusAccessor(self)) def apply_run_settings( self, @@ -235,7 +450,7 @@ def run_status_convenience(Obj): as_dataframe=False, **kwargs, ): - stat = self.run_status( + stat = self._fetch_run_status_dict( run_number=run_number, par_type=par_type, force_reload=force_reload, diff --git a/eco/timing/timing_diag.py b/eco/timing/timing_diag.py index eb844a9..6be563d 100644 --- a/eco/timing/timing_diag.py +++ b/eco/timing/timing_diag.py @@ -23,6 +23,7 @@ from pathlib import Path import datahub as dh from pandas import DataFrame from scipy.optimize import curve_fit +import pickle # from time import sleep @@ -395,7 +396,8 @@ class TimetoolBerninaUSD(Assembly): y, ymed, yerr, - to_display=True, + pl=None, + filepath="", to_elog=True, path_figure="", filepath_data="", @@ -413,10 +415,19 @@ class TimetoolBerninaUSD(Assembly): ax1 = fig.add_subplot(gs[0, 1:4], sharey=ax0) ax2 = fig.add_subplot(gs[0, 4:]) ax1.pcolor(1e15 * x, bins_center, hists) - line = ax1.errorbar(1e15 * x, ymed, yerr, color="red", marker=".", linestyle="") + line = ax1.errorbar( + 1e15 * np.asarray(x), ymed, yerr, color="red", marker=".", linestyle="" + ) fit = ax1.plot( 1e15 * np.polyval(p, ymed), ymed, label="poly fit", color="yellow" ) + if pl is not None: + fitl = ax1.plot( + 1e15 * np.polyval(pl, ymed), + ymed, + label=f"poly fit last calibration\n{filepath.stem}", + color="orange", + ) ax0.axvline(0, linestyle="--", color="k") ax0.plot(1e15 * (np.polyval(p, ymed) - x), ymed, color="royalblue") ax1.set_xlabel("tt_kb.delay (fs)") @@ -433,12 +444,18 @@ class TimetoolBerninaUSD(Assembly): ax2.hist(at, bins="auto", edgecolor="black", histtype="step") try: fx = d[1][:-1] + d[1][1] - d[1][0] - parsopt, parss = curve_fit(self.gauss, fx, d[0], [30, 0, 50]) + parsopt, parss = curve_fit( + self.gauss, + fx, + d[0], + [30, 0, 50], + bounds=[[5, -1000, 5], [500, 1000, 50000]], + ) plx = np.arange(np.min(d[1]), np.max(d[1]), 0.1) ax2.plot(plx, self.gauss(plx, *parsopt), color="k") except: print("Fitting of arrival time histogram failed") - parsopt = [0, 0, 0] + parsopt = [0.0, 0.0, 0.0] pass ax0.set_title("Residual") ax1.set_title("Scan") @@ -485,6 +502,25 @@ class TimetoolBerninaUSD(Assembly): except Exception as e: print(f"Elog posting failed with:\n {e}") + def load_last_calib(self, datapath): + files = sorted(Path(datapath).glob("*_calib.pkl")) + if not files: + raise FileNotFoundError(f"No previous calibration files found in {datapath}") + with open(files[-1], "rb") as file: + lc = pickle.load(file) + return lc["p"], files[-1] + + def save_calibration(self, p, x, y, ymed, yerr, dpath_calib): + lc = { + "p": p, + "tt_kb.delay": x, + "tt_kb.edge_position_px": y, + "ymed": ymed, + "yerr": yerr, + } + with open(dpath_calib, "wb") as file: + pickle.dump(lc, file) + def calibrate( self, seconds=5, @@ -504,14 +540,9 @@ class TimetoolBerninaUSD(Assembly): timestamp = datetime.now().strftime("%Y-%m-%d_%H-%M-%S") pgroup = config_bernina.pgroup() - path_data = ( - "/sf/bernina/config/src/beamline_devices/tt_kb/data/" - + f"{timestamp}_{pgroup}" - ) - path_figure = ( - "/sf/bernina/config/src/beamline_devices/tt_kb/figures/" - + f"{timestamp}_{pgroup}" - ) + basepath = "/sf/bernina/config/src/beamline_devices/tt_kb/" + path_data = f"{basepath}data/{timestamp}_{pgroup}" + path_figure = f"{basepath}figures/{timestamp}_{pgroup}" feedback = self.feedback_enabled() t0 = self.delay() @@ -550,15 +581,65 @@ class TimetoolBerninaUSD(Assembly): x, df, pids_start, pids_stop, filepath=path_data + ".esc.h5" ) + ########## scan backwards ########## + if bidirectional: + xb, pids_start, pids_stop = self.scan_calibration( + seconds=seconds, + scan_range=scan_range, + scan_steps=scan_steps, + reverse_direction=not reverse_direction, + ) + + yb, dfb = self.retrieve_calibration_data( + pids_start=pids_start, + pids_stop=pids_stop, + additional_channels=additional_channels, + ) + self.dataframe_to_escape_dataset( + xb, dfb, pids_start, pids_stop, filepath=path_data + ".esc.h5" + ) + x = np.concatenate([x, xb]) + y = y + yb + ########## fit data ########## p, x, y, ymed, yerr = self.fit_calibration_data( x, y, filter_outliers=filter_outliers ) + ####### load last calib (before saving the new one!) ###### + pl = None + filepath = None + edgel = None + try: + pl, filepath = self.load_last_calib(basepath + "data/") + print(f"Compare new calibration to last calibration taken: \n{filepath}") + for root in np.roots(pl): + if np.isreal(root) and 0 < root.real < 2000: + edgel = root.real + break + except Exception as e: + print("Failed to load last calibration") + print(e) + + # User question: keep pixel of previous calib + if edgel is not None: + ans = "" + while not any([a in ans for a in ["y", "n"]]): + try: + ans = input( + f"Do you wish to shift the calibration to keep the edge at the same pixel ({edgel:.5}) as in the previous calibration (y/n)?" + ) + except: + continue + if ans == "y": + p[-1] = -(p[0] * edgel**2 + p[1] * edgel) + print(f"Shifted calibration curve to preserve edge position: {p}") + elif ans == "n": + continue + ####### save calibration ###### dpath_calib = Path(f"{path_data}_calib.pkl") - df = DataFrame({"tt_kb.delay": x, "tt_kb.edge_position_px": y}) - df.to_pickle(dpath_calib) + self.save_calibration(p, x, y, ymed, yerr, dpath_calib) ########## plot data ########## if plot: @@ -568,6 +649,8 @@ class TimetoolBerninaUSD(Assembly): y, ymed, yerr, + pl, + filepath, to_elog=to_elog, path_figure=path_figure, filepath_data=dpath_calib, @@ -575,6 +658,18 @@ class TimetoolBerninaUSD(Assembly): if update_pipeline_config: self.set_calibration_values(p, pipeline=pipeline, to_elog=to_elog) + else: + # User question: apply calib + ans = "" + while not any([a in ans for a in ["y", "n"]]): + try: + ans = input(f"Do you wish to update the pipeline config (y/n)?") + except: + continue + if ans == "y": + self.set_calibration_values(p, pipeline=pipeline, to_elog=to_elog) + elif ans == "n": + continue if feedback: self.feedback_enabled(1)