timetool calib work and scan correctness

This commit is contained in:
gac bernina
2026-07-11 14:08:43 +02:00
parent a8945f22c1
commit 4899bf0bd8
3 changed files with 387 additions and 41 deletions
+60 -24
View File
@@ -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,
+218 -3
View File
@@ -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: <live value>)]"
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[<TAB>]`. 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[<TAB>]."""
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"<run_status for {self._obj.alias.get_full_name()}: {avail}>\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,
+109 -14
View File
@@ -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)