Files
Jungfraujoch/acquisition_device/AcquisitionDevice.cpp
T
leonarski_fandClaude Opus 5 16bf3408f0 Address code-review findings; make detection limits detector-driven
One changeset, developed together in response to a review of this branch, so the
files carry several of the changes at once. Full test suite passes (733 cases).

Spot finding
- Split ImageSpotFinder into Detect() (flag strong pixels - the expensive
  per-pixel pass) and ExtractSpots() (CCL + min/max-pix + resolution mask), with
  Run() = both. The per-image min-pix escalation now detects ONCE and repeats
  only the cheap extraction, instead of re-running the whole finder four times
  per frame as it did on the default path. It also keeps the winning attempt's
  spot list rather than re-extracting it, so the frame that is integrated is
  exactly the frame that was scored - which a GPU re-extract could not guarantee
  (float atomic ordering).
- spot_finding_time_s no longer swallows indexing time, and indexing_time_s now
  sums every escalation call instead of reporting only the last.

Detection limits follow the detector
- The azimuthal-integration upper q and the spot-finding high-resolution limit
  are now std::optional, in the C++ structs AND in the OpenAPI schema, and
  resolve to the detector's own maximum (DiffractionExperiment::GetDetectorMaxQ_
  recipA). Adaptive detection reads a pixel's ring from the azimuthal bins, so a
  pixel outside that q range could never be strong - the integration range
  silently bounded what detection could see, regardless of the requested
  resolution limit. Regenerated the C++ and TypeScript clients; the viewer and
  the web frontend each gained a "to detector edge" switch.

Detection defaults are now per workflow (measured, not assumed)
- Stills: adaptive detection, min-pix chosen per image, no resolution clipping.
- Rotation: fixed-threshold finder, min-pix 2, 1.5 A limit.
  On a 33-crystal rotation battery, adaptive detection helped four hard crystals
  but deterministically broke three (a lost space group, a halved indexing rate,
  a collapsed merge), and the detector-edge limit cost indexing on a strong
  rotation set (100.0 -> 96.8%). Each is still overridable by its flag, and
  --no-adaptive-spots is new.

Indexer seed escalation
- Stop escalating once a seed's lattice explains >= 90% of the seed spots.
  Previously any frame with >= 80 spots always paid three indexer calls, online
  broker included.

Merge-consistency filter
- --min-image-cc gated on a per-image CC computed BEFORE the stills partiality
  post-refinement and never refreshed; the refiner now recomputes it, so the
  reported CC describes the data that are actually merged.
- Replaced the per-call cc_mask argument with one MergeOnTheFly flag, so the
  merge, the error model and MergeStats can no longer disagree about which
  images are in (the --scale path merged unfiltered while its statistics were
  filtered).

Per-image B-factor refinement (-B) removed
- Measured on four serial-stills datasets: it is a no-op where the per-image fit
  is well conditioned and actively harmful where it is not (CC1/2 -8.1, R_meas
  +23.2 on the weakest large-cell set, whose fits hit their [-50, 200] bounds on
  14-25% of images). It had also been silently DISCARDED since the partiality
  post-refinement landed - reported but not applied. Rather than fix and keep a
  knob with no demonstrated benefit, the flag and the whole image_scale_b_factor
  chain are gone: setting, scaling fit, message field, CBOR, HDF5 write and
  read-back, per-image plot, OpenAPI enum, viewer column and checkbox, docs.
  ScaleOnTheFly no longer needs Ceres at all - the fit is a linear IRLS.
  (The Wilson per-image b_factor is a different quantity and stays.)

Stills partiality width now fits both of its components
- sigma^2 = gamma0^2 + (gamma_e*d*)^2 instead of a purely angular gamma_e*d*
  with gamma0 pinned to 0. Fitted per crystal by least squares of dist_ewald^2
  on d*^2. The angular-only width is fitted over a d*^2-dense population, so it
  was pinned by the high-resolution edge and collapsed at low d*: median
  partiality 0.008 beyond 13 A for reflections that were plainly recorded, 55%
  of them under the merge's partiality floor, and the survivors divided by those
  values - which inflated the merged low-resolution intensity scale 3.6x
  (~ +9 A^2 of apparent B). Measured on 5000 stills: the ramp flattens to 0.89x,
  no observation is dropped any more (701750 -> 716811), shell-mean CC1/2 and
  R-free improve slightly. Note CC1/2, R_meas, completeness and a B-refining
  R-free are all blind to that ramp, which is why it survived earlier validation;
  the cost is high-resolution R_meas (98.5 -> 101.9 shell-averaged).

Removed dead code from add-then-remove churn
- Prediction-time "still partiality" (unreachable: no setter), the phantom
  IndexingSettings::min_indexed_spot_fraction knob (getter, no setter - now the
  constant it always was), StillsPartialityRefine's caller-less Settings
  constructor and its reference to a long-gone env var, ProcessImage's unread
  bool return, an unused include, and a dead viewer overlay hook.

Also
- Viewer: the magnifier compared a QImage with itself, so its scene rect was set
  once ever and it could not pan into a larger dataset; the hover tail timer
  could fire after leaveEvent and resurrect the resolution readout outside the
  image.
- update_version.sh regenerated the frontend lock file BEFORE bumping the
  version (every release shipped an off-by-one lock), and did git rm/git add on
  a path that has not existed since the client moved to src/client - with no
  set -e, both failed silently.
- fpga/pcie_driver/postinstall.sh tested "[ ! occurrences > 0 ]", which is a
  redirect, not a test, so dkms add never ran.
- Unit tests for the adaptive-threshold host functions, which had none.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
2026-07-27 09:07:00 +02:00

330 lines
13 KiB
C++

// SPDX-FileCopyrightText: 2024 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
// SPDX-License-Identifier: GPL-3.0-only
#include <thread>
#include <fstream>
#include <cmath>
#include "../common/JFJochException.h"
#include "AcquisitionDevice.h"
#include "../common/NetworkAddressConvert.h"
AcquisitionDevice::AcquisitionDevice(uint16_t in_data_stream) {
logger = nullptr;
data_stream = in_data_stream;
}
void AcquisitionDevice::PrepareAction(const DiffractionExperiment &experiment) {
if (experiment.GetModulesNum(data_stream) > max_modules)
throw(JFJochException(JFJochExceptionCategory::InputParameterAboveMax,
"Number of modules exceeds max possible for FPGA"));
counters.Reset(experiment, data_stream);
}
void AcquisitionDevice::StartAction(const DiffractionExperiment &experiment, uint32_t optional_flags) {
Cancel();
if (experiment.GetModulesNum(data_stream) > max_modules)
throw(JFJochException(JFJochExceptionCategory::InputParameterAboveMax,
"Number of modules exceeds max possible for FPGA"));
counters.Reset(experiment, data_stream);
expected_frames = experiment.GetFrameNum() / experiment.GetFPGASummation();
// Ensure internal WR queue is empty
work_request_queue.Clear();
Start(experiment, optional_flags);
for (uint32_t i = 0; i < buffer_device.size(); i++)
SendWorkRequest(i);
auto c = work_completion_queue.GetBlocking();
if (c.type != Completion::Type::Start)
throw JFJochException(JFJochExceptionCategory::AcquisitionDeviceError, "Mismatch in work completions");
StartSendingWorkRequests();
start_time = std::chrono::system_clock::now();
if (experiment.IsUsingInternalPacketGen())
RunInternalGenerator(experiment);
}
void AcquisitionDevice::WaitForActionComplete() {
auto c = work_completion_queue.GetBlocking();
while (c.type != Completion::Type::End) {
DeviceOutput* output;
bool skip_analysis = false;
try {
output = GetDeviceOutput(c.handle);
} catch (const JFJochException &e) {
if (logger)
logger->ErrorException(e);
skip_analysis = true;
}
if (!skip_analysis) {
c.module_number = output->module_statistics.module_number;
c.packet_count = output->module_statistics.packet_count;
c.frame_number = output->module_statistics.frame_number;
if (c.frame_number >= expected_frames) {
Cancel();
// this frame is not of any interest, therefore its location can be immediately released
SendWorkRequest(c.handle);
} else if (c.module_number >= max_modules) {
// Module number out of bounds, don't process
if (logger != nullptr)
logger->Error("Completion with wrong module number data stream {} completion frame number {} module {} handle {}",
data_stream, c.frame_number, c.module_number, c.handle);
SendWorkRequest(c.handle);
} else if (c.frame_number < counters.GetSlowestFrameNumber()) {
// Module is falling behind, needs to return the handle then
SendWorkRequest(c.handle);
} else {
try {
counters.UpdateCounters(&c);
} catch (const JFJochException &e) {
if (logger)
logger->ErrorException(e);
SendWorkRequest(c.handle);
}
}
if (logger != nullptr)
logger->Debug("Data stream {} completion frame number {} module {} handle {}",
data_stream, c.frame_number, c.module_number, c.handle);
}
c = work_completion_queue.GetBlocking();
}
counters.SetAcquisitionFinished();
end_time = std::chrono::system_clock::now();
Cancel();
Finalize();
}
void AcquisitionDevice::SendWorkRequest(uint32_t handle) {
work_request_queue.Put(WorkRequest{
.handle = handle
});
}
uint64_t AcquisitionDevice::GetBytesReceived() const {
return counters.GetBytesReceived();
}
const DeviceOutput *AcquisitionDevice::GetDeviceOutput(size_t frame_number, uint16_t module_number) const {
auto handle = counters.GetBufferHandle(frame_number, module_number);
if (handle != HandleNotValid)
return GetDeviceOutput(handle);
else
throw JFJochException(JFJochExceptionCategory::ArrayOutOfBounds, "Frame not collected");
}
const DeviceOutput *AcquisitionDevice::GetDeviceOutput(size_t handle) const {
if (handle >= buffer_device.size())
throw JFJochException(JFJochExceptionCategory::ArrayOutOfBounds, "Handle outside of range");
else
return (DeviceOutput *) buffer_device.at(handle);
}
DeviceOutput *AcquisitionDevice::GetDeviceOutput(size_t handle) {
if (handle >= buffer_device.size())
throw JFJochException(JFJochExceptionCategory::ArrayOutOfBounds, "Handle outside of range");
else
return (DeviceOutput *) buffer_device.at(handle);
}
void AcquisitionDevice::InitializeCalibration(const DiffractionExperiment &experiment, const JFCalibration &calib) {}
void AcquisitionDevice::InitializeIntegrationMap(const DiffractionExperiment &experiment,
const std::vector<uint16_t> &v,
const std::vector<float> &weights) {}
void AcquisitionDevice::InitializeIntegrationMap(const uint16_t *map, const float *weights, size_t module_number) {}
void AcquisitionDevice::InitializeSpotFinderResolutionMap(const float *data, size_t module_number) {}
void AcquisitionDevice::InitializeROIMap(const uint16_t *map, size_t module_number) {}
void AcquisitionDevice::InitializePixelMask(const uint32_t *module_mask, size_t module_number) {}
void AcquisitionDevice::InitializeROIMap(const DiffractionExperiment& experiment, const std::vector<uint16_t>& roi_map) {
if (roi_map.size() != experiment.GetXPixelsNumConv() * experiment.GetYPixelsNumConv())
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "Mismatch in array size");
std::vector<uint16_t> tmp(RAW_MODULE_SIZE);
auto offset = experiment.GetFirstModuleOfDataStream(data_stream);
size_t modules = experiment.GetModulesNum(data_stream);
for (int m = 0; m < modules; m++) {
ConvertedToRawGeometry(experiment, offset + m, tmp.data(), roi_map.data());
InitializeROIMap(tmp.data(), m);
}
}
void AcquisitionDevice::InitializeEmptyPixelMask(const DiffractionExperiment &experiment) {
std::vector<uint32_t> empty_mask(RAW_MODULE_SIZE);
size_t modules = experiment.GetModulesNum(data_stream);
for (int m = 0; m < modules; m++)
InitializePixelMask(empty_mask.data(), m);
}
void AcquisitionDevice::InitializeDataProcessing(const DiffractionExperiment &experiment,
const AzimuthalIntegrationMapping &azint) {
auto offset = experiment.GetFirstModuleOfDataStream(data_stream);
size_t modules = experiment.GetModulesNum(data_stream);
// When azimuthal integration is forced onto the CPU, the FPGA must not bin pixels
// (the integration map can address more bins than the FPGA supports); the CPU path
// computes the profile from the assembled image instead.
const bool load_integration_map = !experiment.GetAzimuthalIntegrationSettings().IsForceCPUinFPGAWorkflow();
if (experiment.IsGeometryTransformed()) {
std::vector<float> tmp1(RAW_MODULE_SIZE);
std::vector<uint16_t> tmp2(RAW_MODULE_SIZE);
for (int m = 0; m < modules; m++) {
if (load_integration_map) {
ConvertedToRawGeometry(experiment, offset + m, tmp1.data(), azint.Corrections().data());
ConvertedToRawGeometry(experiment, offset + m, tmp2.data(), azint.GetPixelToBin().data());
InitializeIntegrationMap(tmp2.data(), tmp1.data(), m);
}
ConvertedToRawGeometry(experiment, offset + m, tmp1.data(), azint.Resolution().data());
InitializeSpotFinderResolutionMap(tmp1.data(), m);
}
} else {
for (int m = 0; m < modules; m++) {
if (load_integration_map)
InitializeIntegrationMap(azint.GetPixelToBin().data() + (offset + m) * RAW_MODULE_SIZE,
azint.Corrections().data() + (offset + m) * RAW_MODULE_SIZE,
m);
InitializeSpotFinderResolutionMap(azint.Resolution().data() + (m + offset) * RAW_MODULE_SIZE,
m);
}
}
}
void AcquisitionDevice::InitializePixelMask(const DiffractionExperiment &experiment, const PixelMask &mask) {
auto offset = experiment.GetFirstModuleOfDataStream(data_stream);
size_t modules = experiment.GetModulesNum(data_stream);
std::vector<uint32_t> tmp(RAW_MODULE_SIZE);
for (int m = 0; m < modules; m++) {
ConvertedToRawGeometry(experiment, offset + m, tmp.data(), mask.GetMask().data());
InitializePixelMask(tmp.data(), m);
}
}
void AcquisitionDevice::FrameBufferRelease(size_t frame_number, uint16_t module_number) {
auto handle = counters.GetBufferHandleAndClear(frame_number, module_number);
if (handle != AcquisitionCounters::HandleNotFound)
SendWorkRequest(handle);
}
void AcquisitionDevice::EnableLogging(Logger *in_logger) {
logger = in_logger;
}
int32_t AcquisitionDevice::GetNUMANode() const {
return -1;
}
uint16_t AcquisitionDevice::GetUDPPort() const {
return 1234;
}
const AcquisitionCounters &AcquisitionDevice::Counters() const {
return counters;
}
std::string AcquisitionDevice::GetIPv4Address() const {
return IPv4AddressToStr(ipv4_addr);
}
std::string AcquisitionDevice::GetMACAddress() const {
return MacAddressToStr(mac_addr);
}
DataCollectionStatus AcquisitionDevice::GetDataCollectionStatus() const {
return {};
}
DeviceStatus AcquisitionDevice::GetDeviceStatus() const {
return {};
}
AcquisitionDeviceStatistics AcquisitionDevice::GetStatistics() const {
AcquisitionDeviceStatistics ret{};
ret.bytes_received = GetBytesReceived();
ret.start_timestamp = std::chrono::system_clock::to_time_t(start_time);
ret.end_timestamp = std::chrono::system_clock::to_time_t(end_time);
ret.packets_expected = counters.GetTotalExpectedPackets();
ret.good_packets = counters.GetTotalPackets();
for (int i = 0; i < counters.GetModuleNumber(); i++)
ret.packets_received_per_module.push_back(counters.GetTotalPackets(i));
if ((ret.packets_expected == 0) || (ret.good_packets == ret.packets_expected))
ret.efficiency = 1.0;
else
ret.efficiency = static_cast<float>(ret.good_packets) / static_cast<float>(ret.packets_expected);
return ret;
}
void AcquisitionDevice::SetIPv4Address(uint32_t ipv4_addr_network_order) {
ipv4_addr = ipv4_addr_network_order;
}
AcquisitionDeviceNetConfig AcquisitionDevice::GetNetConfig() const {
return {
.mac_addr = GetMACAddress(),
.ipv4_addr = GetIPv4Address(),
.udp_port = GetUDPPort()
};
}
void AcquisitionDevice::RunInternalGenerator(const DiffractionExperiment &experiment) {
FrameGeneratorConfig config{};
config.frames = experiment.GetFrameNum() + DELAY_FRAMES_STOP_AND_QUIT + 1;
config.modules = experiment.GetModulesNum(data_stream);
config.data_stream = data_stream;
config.pulse_id = INT_PKT_GEN_BUNCHID;
config.exptime = INT_PKT_GEN_EXPTTIME;
config.debug = INT_PKT_GEN_DEBUG;
config.dest_mac_addr = MacAddressFromStr(GetMACAddress());
config.dest_ipv4_addr = IPv4AddressFromStr(GetIPv4Address());
config.images_in_memory = experiment.GetInternalPacketGeneratorImages() - 1;
switch (experiment.GetDetectorSetup().GetDetectorType()) {
case DetectorType::JUNGFRAU:
config.detector_type = SLS_DETECTOR_TYPE_JUNGFRAU;
break;
case DetectorType::EIGER:
config.detector_type = SLS_DETECTOR_TYPE_EIGER;
config.eiger_bit_depth = experiment.GetBitDepthReadout();
break;
default:
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "Detector not supported");
}
HW_RunInternalGenerator(config);
}
void AcquisitionDevice::SetSpotFinderParameters(const SpotFindingSettings &settings) {
SpotFinderParameters fpga_parameters{};
fpga_parameters.snr_threshold = settings.signal_to_noise_threshold;
fpga_parameters.count_threshold = settings.photon_count_threshold;
fpga_parameters.max_d = settings.low_resolution_limit;
fpga_parameters.min_d = settings.high_resolution_limit.value_or(0.0f);
fpga_parameters.min_pix_per_spot = settings.min_pix_per_spot.value_or(2);
HW_SetSpotFinderParameters(fpga_parameters);
}