Merge branch 'report-work' into rc175: beam-centre warning, overloads and overlaps, XDS geometry keywords in the report
Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01SVmAWnzCmRKAXVUCdc4iNi
This commit is contained in:
leonarski_fandClaude Opus 5.5 committed
commit
bc9f2caccb
14 files changed
+509
-88
No files matched your search
+14
-10
@@ -504,23 +504,27 @@ sweep that does not have them will get this far and no further.
|
||||
|
||||
## Comparing the geometry with XDS
|
||||
|
||||
Every run logs the detector geometry a second time in **XDS's convention**, so it can be read
|
||||
straight across against the `IDXREF.LP` / `CORRECT.LP` of an XDS run on the same data:
|
||||
The report states the geometry the result was integrated at in **XDS's convention**, in its
|
||||
`XDS GEOMETRY KEYWORDS` section (see [the results report](RUGNUX_REPORT.md)), so it can be read
|
||||
straight across against the `IDXREF.LP` / `CORRECT.LP` of an XDS run on the same data, or pasted into
|
||||
`XDS.INP`:
|
||||
|
||||
```
|
||||
XDS convention: ORGX= 1091.00 ORGY= 1137.00 DETECTOR_DISTANCE= 75.0000
|
||||
XDS convention: DIRECTION_OF_DETECTOR_X-AXIS= 1.000000 0.000000 0.000000
|
||||
XDS convention: DIRECTION_OF_DETECTOR_Y-AXIS= 0.000000 1.000000 0.000000
|
||||
XDS convention: INCIDENT_BEAM_DIRECTION= 0 0 1 X-RAY_WAVELENGTH= 1.000000 QX= QY= 0.075000
|
||||
XDS convention: ROTATION_AXIS= -1.000000 0.000000 0.000000
|
||||
# XDS.INP> ORGX= 1091.00 ORGY= 1137.00
|
||||
# XDS.INP> DETECTOR_DISTANCE= 75.000
|
||||
# XDS.INP> INCIDENT_BEAM_DIRECTION= 0 0 1 X-RAY_WAVELENGTH= 1.000000
|
||||
# XDS.INP> ROTATION_AXIS= 1.000000 0.000000 0.000000
|
||||
```
|
||||
|
||||
XDS is never given this geometry — the [XDS plugin](SOFTWARE_INTEGRATION.md) supplies image data
|
||||
only, and XDS refines its own from `XDS.INP` — which is what makes the comparison worth having. The
|
||||
The [XDS plugin](SOFTWARE_INTEGRATION.md) supplies image data only, so XDS refines its own geometry
|
||||
from `XDS.INP` - which is what makes the comparison worth having. The
|
||||
two laboratory frames coincide (x along increasing detector column, y along increasing row, z along
|
||||
the beam), so the numbers are directly comparable, and a tilt appears as the two detector axis
|
||||
vectors rather than as angles, which is how XDS reports it after refinement. Two things to keep in
|
||||
mind: **`ORGX`/`ORGY` are 1-based**, because XDS counts pixels from 1 and Jungfraujoch from 0; and
|
||||
they are the **PONI**, the same quantity Jungfraujoch's beam centre is — so no correction is needed —
|
||||
but not the direct beam once the detector is tilted (see above).
|
||||
but not the direct beam once the detector is tilted (see above). **`ROTATION_AXIS` is rugnux's axis
|
||||
negated**: rugnux turns the crystal by -phi about its axis, XDS by +phi, and with rugnux's own sign
|
||||
`IDXREF` indexes only the chance floor of a mis-signed spindle. The detector axes appear only where the
|
||||
detector is tilted.
|
||||
|
||||
+87
-2
@@ -128,6 +128,57 @@ JFJOCH_DATASET_SETTINGS= {"beam_x_pxl": 2078.24, "beam_y_pxl": 2233.92, "detecto
|
||||
grep '^JFJOCH_DATASET_SETTINGS=' out_report.txt | cut -d' ' -f2- > geometry.json
|
||||
```
|
||||
|
||||
**The beam centre in the file.** Section 1 sets the file's centre beside the two the run measured:
|
||||
**`FILE_DIRECT_BEAM=`** is where the input file puts the direct beam (its beam centre, under its own
|
||||
detector tilt), **`BACKGROUND_DIRECT_BEAM=`** where the scattered background puts it, with the fit's
|
||||
sigma (the pre-scan's beam-centre check), and **`BEAM_CENTRE_SHIFT=`** how far the centre the result was
|
||||
integrated at (`DIRECT_BEAM=`, the refined one on a rotation run) lies from the file's, in pixels.
|
||||
**`BEAM_CENTRE_TOLERANCE=`** is how far this crystal on this geometry lets the centre be wrong - the
|
||||
displacement at which the accumulated reflections of the longest cell axis smear to 70% of their
|
||||
height, the number the log's `Beam centre check` line states, here on the run's own cell. Where the
|
||||
shift is more than three times that, a **`BEAM_CENTRE_OFF`** warning gives the file's, the refined and the background
|
||||
centre, the difference in pixels and millimetres, and the value to correct the beamline (or the
|
||||
header) to: the run itself integrated at the refined centre, but every other program given the same
|
||||
file starts from the wrong one. A run that refined no centre (stills) is warned on the background
|
||||
alone, where its fit resolves the difference (three sigma) and the difference exceeds three times the
|
||||
tolerance, worded as a check. Between one and three times the tolerance the summary's `Beam centre`
|
||||
line says so and recommends correcting the beamline value, without a warning - most headers are a few
|
||||
pixels off, and at one times the tolerance the warning fired on 107 of 250 battery sets.
|
||||
|
||||
```
|
||||
FILE_DIRECT_BEAM= 1564.31 1702.91
|
||||
BACKGROUND_DIRECT_BEAM= 1571.30 1699.85 0.62
|
||||
BEAM_CENTRE_SHIFT= 7.15
|
||||
BEAM_CENTRE_TOLERANCE= 3.49
|
||||
```
|
||||
|
||||
**`XDS GEOMETRY KEYWORDS`**, the section after section 1, is the same geometry as the `XDS.INP`
|
||||
keywords that carry it, ready to paste into an XDS run: `ORGX`/`ORGY` (the PONI, numbered from 1 as XDS
|
||||
numbers pixels), `DETECTOR_DISTANCE`, the detector axes `DIRECTION_OF_DETECTOR_X-AXIS`/`Y-AXIS` where the
|
||||
detector is tilted (they carry the tilt), `INCIDENT_BEAM_DIRECTION` with `X-RAY_WAVELENGTH`, and
|
||||
`ROTATION_AXIS` - the axis the run integrated at, **negated**, because XDS turns the crystal the other way
|
||||
about the axis it is given (checked with XDS `IDXREF`). Geometry only: XDS determines the cell and the
|
||||
space group itself. The lines carry a fixed prefix, so
|
||||
|
||||
```bash
|
||||
sed -n 's/^# XDS.INP> //p' myrun_report.txt
|
||||
```
|
||||
|
||||
prints them as XDS reads them:
|
||||
|
||||
```
|
||||
ORGX= 749.15 ORGY= 817.32
|
||||
DETECTOR_DISTANCE= 175.668
|
||||
DIRECTION_OF_DETECTOR_X-AXIS= 1.000000 0.000000 0.000444
|
||||
DIRECTION_OF_DETECTOR_Y-AXIS= 0.000000 1.000000 -0.000976
|
||||
INCIDENT_BEAM_DIRECTION= 0 0 1 X-RAY_WAVELENGTH= 0.918400
|
||||
ROTATION_AXIS= 1.000000 0.000000 0.000000
|
||||
```
|
||||
|
||||
The conventions are XDS's own (W. Kabsch, "XDS" (2010), Acta Cryst. D66, 125-132
|
||||
[doi:10.1107/S0907444909047337](https://doi.org/10.1107/S0907444909047337)); this section is the one
|
||||
place the run states them.
|
||||
|
||||
**Which pass.** A rotation run integrates twice — once at the geometry in the input file (the first
|
||||
half of the sweep but at least 120°, `--prepass-fraction`), then the whole sweep again at the post-refined geometry — and can integrate a third time if a guard rejects the second pass.
|
||||
There is **one** report, for the pass that became the canonical output, and `PASS=` /
|
||||
@@ -156,13 +207,13 @@ and a run that produced a textbook data set read identically for their first thr
|
||||
`GONIO_SCALE`, `SPINDLE_CAP`, `ANISOTROPY`, `TWINNING`, `PSEUDO_TRANSLATION`,
|
||||
`LATTICE_TRANSLATION`, `MODEL_HAND`, `MODEL_NOT_VALIDATED`, `REFERENCE_MISMATCH`, `CANCELLED`, `RESOLUTION_FIT`,
|
||||
`FLIGHT_PATH`, `GEOMETRY_NOT_CONVERGED`, `SCALING_NOT_CONVERGED`, `HARMONIC_CONTAMINATION`,
|
||||
`SUPERCELL_POSSIBLE`, `MULTIPLE_LATTICES`, `ICE_RINGS`, `POWDER_RINGS`. `NONE` when nothing fired. A code appears if and only if its
|
||||
`SUPERCELL_POSSIBLE`, `MULTIPLE_LATTICES`, `ICE_RINGS`, `POWDER_RINGS`, `BEAM_CENTRE_OFF`, `OVERLOADS`, `OVERLAPS`. `NONE` when nothing fired. A code appears if and only if its
|
||||
warning fired, so the flags and the `WARNING:` lines are two renderings of one list. Thresholds are
|
||||
deliberately low: a warning is a prompt to check something, not a verdict, and some of them fire on
|
||||
data that turn out fine — the closed
|
||||
type for machinery, the open sentence for a person.
|
||||
- Below them, one line each, the facts a reader needs before reading further: space group, cell,
|
||||
mosaicity, the powder and ice rings, multiple lattices, resolution, completeness, signal, anomalous
|
||||
mosaicity, the beam centre, the powder and ice rings, overloads, overlapping spots, multiple lattices, resolution, completeness, signal, anomalous
|
||||
signal, anisotropy, pseudo-symmetry, twinning, supercell, radiation damage and the sweep. A line
|
||||
states a finding that is not a condition too - weak powder rings, a possible supercell no stronger than
|
||||
correct cells show, an arbitrary indexing choice (`Indexing choice`, a rotation sweep whose cell
|
||||
@@ -497,6 +548,40 @@ scaling. The summary has an `Ice` line, and where the merge found ice an `ICE_RI
|
||||
look at the ice-ring shells (one warning line, whichever of the two measurements raised it). The measurement is described in
|
||||
[CPU/GPU data analysis ▸ Resolution and ice-ring handling](CPU_DATA_ANALYSIS_IMAGE.md#33-resolution-and-ice-ring-handling).
|
||||
|
||||
## Overloads and overlapping spots
|
||||
|
||||
Two ways the measurement, rather than the crystal, cost reflections - each with what to change at the
|
||||
next collection. Both are in section 4, both have a summary line, and both warn at a low threshold.
|
||||
|
||||
```
|
||||
OVERLOADED_REFLECTIONS= 1740
|
||||
OVERLOADED_REFLECTIONS_PCT= 0.017
|
||||
OVERLOADED_REFLECTIONS_LOW_RES_PCT= 0.418
|
||||
OVERLAPPED_REFLECTIONS_PCT= 11.10
|
||||
OVERLAPPED_REFLECTIONS_MAJOR_PCT= 0.44
|
||||
```
|
||||
|
||||
(a strongly diffracting lysozyme sweep: an `OVERLOADS` warning, no `OVERLAPS` one)
|
||||
|
||||
**`OVERLOADED_REFLECTIONS=`** counts the integrated reflections - every frame's partial, before the
|
||||
merge - with a saturated pixel in the signal disk (or one the detector could not count on that frame),
|
||||
and **`OVERLOADED_REFLECTIONS_PCT=`** their share. They are not measurements: the merge leaves them out
|
||||
(a rotation merge drops the whole rocking event, `OBSERVATIONS_REJECTED_OVERLOAD=` in section 3), and
|
||||
they are the strongest reflections the crystal gave, so they concentrate at low resolution;
|
||||
**`OVERLOADED_REFLECTIONS_LOW_RES_PCT=`** is their share among the reflections coarser than 4 Å. An
|
||||
**`OVERLOADS`** warning fires from 0.05% of all integrated reflections or 0.25% of those coarser than 4 Å,
|
||||
and recommends attenuating the beam or shortening the exposure (or a second, attenuated low-resolution
|
||||
sweep).
|
||||
|
||||
**`OVERLAPPED_REFLECTIONS_PCT=`** is the share of the predicted reflections whose signal disk a nearer
|
||||
neighbour's prediction takes pixels from, as the integrator's overlap treatment (`--overlap`, default
|
||||
`exclude`) assigns them, and **`OVERLAPPED_REFLECTIONS_MAJOR_PCT=`** the share left with less than
|
||||
`--overlap-minpk` (default 75%) of the disk its own - the ones XDS's MINPK would reject. Under
|
||||
`exclude` the shared pixels are left out of the profile fit, which keeps the reflection at a cost in
|
||||
precision. An **`OVERLAPS`** warning fires from 1% in the second number and recommends a longer
|
||||
detector distance, finer rotation slicing or a smaller beam. Neither key is written with
|
||||
`--overlap off`, which assigns no pixel to anyone.
|
||||
|
||||
## Index-2 superstructure
|
||||
|
||||
On rotation data, section 4 reports whether the lattice the run adopted has intensity at half-integer
|
||||
|
||||
@@ -201,6 +201,8 @@ void MXAnalysisWithoutFPGA::Analyze(DataMessage &output,
|
||||
probe_counts.bkg_starved += after.bkg_starved - before.bkg_starved;
|
||||
probe_counts.bkg_starved_by_neighbour += after.bkg_starved_by_neighbour - before.bkg_starved_by_neighbour;
|
||||
probe_counts.profile_fallback += after.profile_fallback - before.profile_fallback;
|
||||
probe_counts.overlapped += after.overlapped - before.overlapped;
|
||||
probe_counts.overlapped_major += after.overlapped_major - before.overlapped_major;
|
||||
return ret;
|
||||
};
|
||||
|
||||
@@ -439,6 +441,8 @@ BraggIntegrationCounts MXAnalysisWithoutFPGA::BraggCounts() const {
|
||||
c.bkg_starved -= probe_counts.bkg_starved;
|
||||
c.bkg_starved_by_neighbour -= probe_counts.bkg_starved_by_neighbour;
|
||||
c.profile_fallback -= probe_counts.profile_fallback;
|
||||
c.overlapped -= probe_counts.overlapped;
|
||||
c.overlapped_major -= probe_counts.overlapped_major;
|
||||
return c;
|
||||
}
|
||||
|
||||
|
||||
@@ -134,6 +134,12 @@ struct BraggIntegrationCounts {
|
||||
// costs data adds 4.1 percentage points on top of a floor of 0.02%.
|
||||
uint64_t bkg_starved_by_neighbour = 0;
|
||||
uint64_t profile_fallback = 0; // the profile fit disagreed with the box sum; the box sum was kept
|
||||
// Signal disks a nearer neighbour's prediction takes pixels from (the overlap ownership, so counted
|
||||
// only while an overlap treatment is on), and of those the ones left with less than
|
||||
// overlap_min_peak of the disk their own - what --overlap reject drops, and XDS's MINPK. Reported
|
||||
// only (the crowding of the pattern); nothing reads them.
|
||||
uint64_t overlapped = 0;
|
||||
uint64_t overlapped_major = 0;
|
||||
};
|
||||
|
||||
// One reflection's extracted intensity, produced by the derived engine and turned into a
|
||||
|
||||
@@ -269,6 +269,11 @@ std::vector<Reflection> BraggIntegrationEngineCPU::RunImpl(const Sampler &img,
|
||||
}
|
||||
}
|
||||
|
||||
if (use_owner && n_own < n_disk) {
|
||||
++counts.overlapped;
|
||||
if (n_own < overlap_min_peak * n_disk) ++counts.overlapped_major;
|
||||
}
|
||||
|
||||
int n_bkg_used = n_bkg; // pixels behind the FINAL background value (trim/clip shrink it)
|
||||
// A masked, untrusted, gapped or overloaded pixel inside the signal disk used to discard the
|
||||
// reflection outright. A profile fit does not need it to: the fit is the amplitude of a
|
||||
|
||||
@@ -49,7 +49,9 @@ constexpr int MOM_STRIDE = 3;
|
||||
constexpr int COUNT_BKG_STARVED = 0;
|
||||
constexpr int COUNT_BKG_STARVED_NEIGHBOUR = 1;
|
||||
constexpr int COUNT_PROFILE_FALLBACK = 2;
|
||||
constexpr int COUNT_SLOTS = 3;
|
||||
constexpr int COUNT_OVERLAPPED = 3;
|
||||
constexpr int COUNT_OVERLAPPED_MAJOR = 4;
|
||||
constexpr int COUNT_SLOTS = 5;
|
||||
|
||||
// --- Mark the r2 signal region of every predicted reflection (race-free: one launch writes one
|
||||
// value). The region is the INNER stencil ellipse in the neighbour's own frame; see the CPU engine.
|
||||
@@ -277,6 +279,10 @@ __global__ void boxsum(const float *px_x, const float *px_y, const float *dd,
|
||||
atomicAdd(counts + COUNT_BKG_STARVED, 1ull);
|
||||
if ((s_full || p.partial_ok) && s_nclean <= 5 && s_nring > 5)
|
||||
atomicAdd(counts + COUNT_BKG_STARVED_NEIGHBOUR, 1ull);
|
||||
if (p.overlap && s_nown < s_ndisk) {
|
||||
atomicAdd(counts + COUNT_OVERLAPPED, 1ull);
|
||||
if (s_nown < p.minpk * s_ndisk) atomicAdd(counts + COUNT_OVERLAPPED_MAJOR, 1ull);
|
||||
}
|
||||
s_bkg = s_accept ? ((double) (long long) s_bkgsum / (double) s_nbkg) : 0.0;
|
||||
s_thr = s_bkg + (double) p.bkg_clip_nsigma * sqrt(fmax(s_bkg, 1.0));
|
||||
// The pixel is an integer, so comparing it against the floor of the threshold accepts
|
||||
@@ -1083,11 +1089,13 @@ std::vector<Reflection> BraggIntegrationEngineGPU::Run(const ImagePreprocessorBu
|
||||
|
||||
BraggIntegrationCounts BraggIntegrationEngineGPU::Counts() const {
|
||||
BraggIntegrationCounts out = counts; // `predicted` is counted on the host, as on the CPU engine
|
||||
unsigned long long host[COUNT_SLOTS] = {0, 0, 0};
|
||||
unsigned long long host[COUNT_SLOTS] = {};
|
||||
cuda_err(cudaMemcpyAsync(host, d_counts, sizeof(host), cudaMemcpyDeviceToHost, *stream));
|
||||
cuda_err(cudaStreamSynchronize(*stream));
|
||||
out.bkg_starved = host[COUNT_BKG_STARVED];
|
||||
out.bkg_starved_by_neighbour = host[COUNT_BKG_STARVED_NEIGHBOUR];
|
||||
out.profile_fallback = host[COUNT_PROFILE_FALLBACK];
|
||||
out.overlapped = host[COUNT_OVERLAPPED];
|
||||
out.overlapped_major = host[COUNT_OVERLAPPED_MAJOR];
|
||||
return out;
|
||||
}
|
||||
@@ -51,6 +51,9 @@ namespace PathologyCode {
|
||||
constexpr const char *SCALING_NOT_CONVERGED = "SCALING_NOT_CONVERGED"; // the per-frame scales hit their cap still moving
|
||||
constexpr const char *POWDER_RINGS = "POWDER_RINGS"; // non-ice powder rings hold many spots
|
||||
constexpr const char *ICE_RINGS = "ICE_RINGS"; // hexagonal-ice rings hold many spots
|
||||
constexpr const char *BEAM_CENTRE_OFF = "BEAM_CENTRE_OFF"; // the file's beam centre is wrong for this geometry
|
||||
constexpr const char *OVERLOADS = "OVERLOADS"; // saturated pixels took reflections
|
||||
constexpr const char *OVERLAPS = "OVERLAPS"; // neighbouring spots share signal pixels
|
||||
}
|
||||
|
||||
// A closed type plus an open description. The type is what machinery switches on; the text is what a
|
||||
|
||||
@@ -42,6 +42,9 @@ namespace {
|
||||
// UNDECLARED_LATTICE_TRANSLATION_PCT and INDEXING_AMBIGUITY_OPERATORS; INDEXING_RATE is a
|
||||
// --developer key on rotation runs; POWDER_ICE_* / POWDER_NON_ICE_SPOT_FRACTION;
|
||||
// the ICE_RINGS and POWDER_RINGS flags.
|
||||
// 16: FILE_DIRECT_BEAM / BACKGROUND_DIRECT_BEAM / BEAM_CENTRE_SHIFT / BEAM_CENTRE_TOLERANCE in
|
||||
// section 1 and the BEAM_CENTRE_OFF flag - the file's beam centre against the refined one;
|
||||
// the XDS GEOMETRY KEYWORDS section; OVERLOADED_* / OVERLAPPED_* and the OVERLOADS / OVERLAPS flags.
|
||||
// 8: every line that is not `KEY= value` and not blank starts with `#`, warnings included
|
||||
// (`# WARNING:`); the header gains RUGNUX_DOWNLOAD and BUILD_CXX_FLAGS, and RUGNUX_GIT is
|
||||
// the full build-time hash.
|
||||
@@ -98,6 +101,22 @@ namespace {
|
||||
// fired on were below it (0.4-0.6% of the sweep).
|
||||
constexpr double SWEEP_GAPS_WARNING_FRACTION = 0.01;
|
||||
|
||||
// The beam centre is warned about beyond this many times what the geometry needs of it
|
||||
// (BEAM_CENTRE_TOLERANCE); between one and this many times the summary states it without a warning.
|
||||
// At one times the tolerance 107 of the 250 battery sets fired - most headers are a few pixels off.
|
||||
constexpr double BEAM_CENTRE_WARNING_NEEDS = 3.0;
|
||||
|
||||
// Overloads are warned about from this share of the integrated reflections, or of those coarser than
|
||||
// OVERLOAD_LOW_RES_A, where a handful already takes the strongest low-resolution terms. Low on purpose:
|
||||
// the strong-diffraction in-house sweep reads 0.017% and 0.42%, and its merge drops 0.12% of its
|
||||
// rocking events; an ordinary sweep with a few hundred overloaded partials reads 0.001% and 0.012%.
|
||||
constexpr double OVERLOAD_WARNING_PCT = 0.05;
|
||||
constexpr double OVERLOAD_LOW_RES_WARNING_PCT = 0.25;
|
||||
|
||||
// Overlapping spots are warned about from this share of the predicted reflections left with less
|
||||
// than --overlap-minpk of their profile their own - the ones XDS's MINPK would reject.
|
||||
constexpr double OVERLAP_WARNING_PCT = 1.0;
|
||||
|
||||
// Whether <|L|> - on the merge reported, or on the one the space-group search was given - reads
|
||||
// outside its physical range. Then no twin verdict is read from it, in either direction.
|
||||
bool TwinStatisticsAbnormal(const ProcessResult &r) {
|
||||
@@ -408,6 +427,8 @@ ReportDocument BuildReportDocument(const std::string &output_prefix,
|
||||
const bool rotation = experiment.IsRotationIndexing();
|
||||
const bool merged = result.has_merge_statistics;
|
||||
doc.warnings = result.warnings;
|
||||
// SUMMARY rows filled by the sections that measure them; empty leaves the row out.
|
||||
std::string beam_centre_row, overload_row, overlap_row;
|
||||
|
||||
// ------------------------------------------------------------------ header block
|
||||
{
|
||||
@@ -543,6 +564,72 @@ ReportDocument BuildReportDocument(const std::string &output_prefix,
|
||||
"integrated at is the best one reached, not a settled one (see PASS_DECISION)");
|
||||
}
|
||||
|
||||
// The beam centre against the file's. The file's is read off the experiment the run was
|
||||
// given - the header, with the header's tilt - and compared with where the result was
|
||||
// integrated (the refined geometry) and with where the scattered background put the beam. The
|
||||
// run corrects a wrong centre for itself; every other program reading the same file does not,
|
||||
// which is why this is worth a warning however well the run went. The tolerance is the one the
|
||||
// pre-scan's check prints (BeamCenterNeed_pxl), taken on the run's own cell.
|
||||
{
|
||||
const auto file_direct = experiment.GetDiffractionGeometry().GetDirectBeam_pxl();
|
||||
const double pixel_mm = experiment.GetPixelSize_mm();
|
||||
const double need = result.beam_centre_need_pxl;
|
||||
const double moved = std::hypot(result.used_direct_beam_x_pxl - file_direct.first,
|
||||
result.used_direct_beam_y_pxl - file_direct.second);
|
||||
const auto &bg = result.background_beam_centre;
|
||||
const double bg_moved = bg ? std::hypot(bg->beam_x_pxl - file_direct.first,
|
||||
bg->beam_y_pxl - file_direct.second) : NAN;
|
||||
Add(s, KeyText("FILE_DIRECT_BEAM", fmt::format("{:.2f} {:.2f}", file_direct.first,
|
||||
file_direct.second)));
|
||||
if (bg)
|
||||
Add(s, KeyText("BACKGROUND_DIRECT_BEAM", fmt::format("{:.2f} {:.2f} {:.2f}", bg->beam_x_pxl,
|
||||
bg->beam_y_pxl, bg->sigma_pxl)));
|
||||
Add(s, KeyReal("BEAM_CENTRE_SHIFT", moved, "{:.2f}"));
|
||||
if (std::isfinite(need))
|
||||
Add(s, KeyReal("BEAM_CENTRE_TOLERANCE", need, "{:.2f}"));
|
||||
// The refined geometry is the evidence wherever the run refined one: the background fit can
|
||||
// be pulled by an asymmetric background, and where both exist and disagree the spot
|
||||
// positions are what the result rests on. Only a run that refined nothing (stills, no
|
||||
// post-refinement) is warned on the background alone, and only where the fit resolves the
|
||||
// difference at all (three of its own sigma, as the pre-scan's check reads it).
|
||||
const bool bg_resolves = bg && bg_moved > 3.0 * bg->sigma_pxl;
|
||||
const std::string bg_text = bg ? fmt::format("the scattered background at ({:.2f}, {:.2f}) +- {:.2f} px",
|
||||
bg->beam_x_pxl, bg->beam_y_pxl, bg->sigma_pxl)
|
||||
: std::string();
|
||||
const double warn_at = BEAM_CENTRE_WARNING_NEEDS * need;
|
||||
if (std::isfinite(need) && moved > warn_at) {
|
||||
beam_centre_row = fmt::format("the file's is OFF by {:.1f} px ({:.2f} mm) - correct the "
|
||||
"beamline value", moved, moved * pixel_mm);
|
||||
Warn(doc, PathologyCode::BEAM_CENTRE_OFF, fmt::format(
|
||||
"the beam centre in the input file is {:.1f} px ({:.2f} mm) from where this run found "
|
||||
"the direct beam, more than three times the {:.1f} px this geometry needs it to be right "
|
||||
"to: the file puts the direct beam at ({:.2f}, {:.2f}), the refined geometry at ({:.2f}, "
|
||||
"{:.2f}){}. This run used the refined centre, but any other program given this file "
|
||||
"starts from the wrong one - correct the beam centre at the beamline (or in the "
|
||||
"file's header) to ({:.2f}, {:.2f})",
|
||||
moved, moved * pixel_mm, need, file_direct.first, file_direct.second,
|
||||
result.used_direct_beam_x_pxl, result.used_direct_beam_y_pxl,
|
||||
bg ? ", " + bg_text : std::string(),
|
||||
result.used_direct_beam_x_pxl, result.used_direct_beam_y_pxl));
|
||||
} else if (std::isfinite(need) && !result.post_refine && bg_resolves && bg_moved > warn_at) {
|
||||
beam_centre_row = fmt::format("the background puts it {:.1f} px ({:.2f} mm) from the file's - "
|
||||
"check the beamline value", bg_moved, bg_moved * pixel_mm);
|
||||
Warn(doc, PathologyCode::BEAM_CENTRE_OFF, fmt::format(
|
||||
"the beam centre in the input file may be off: the file puts the direct beam at "
|
||||
"({:.2f}, {:.2f}), {} - {:.1f} px ({:.2f} mm) away, more than three times the {:.1f} px "
|
||||
"this geometry needs the centre to be right to. Nothing on this run refined the centre, so "
|
||||
"this is the background's word alone: check the beam centre at the beamline and "
|
||||
"correct it if this is confirmed",
|
||||
file_direct.first, file_direct.second, bg_text, bg_moved, bg_moved * pixel_mm, need));
|
||||
} else if (std::isfinite(need) && moved > need)
|
||||
beam_centre_row = fmt::format("the file's is {:.1f} px ({:.2f} mm) from the refined one, more than "
|
||||
"the {:.1f} px this geometry needs but under three times it - worth "
|
||||
"correcting at the beamline", moved, moved * pixel_mm, need);
|
||||
else if (std::isfinite(need))
|
||||
beam_centre_row = fmt::format("the file's is within the {:.1f} px this geometry needs ({:.1f} px "
|
||||
"from the refined one)", need, moved);
|
||||
}
|
||||
|
||||
Add(s, Blank());
|
||||
// Where the geometry came from. This is the whole provenance story a user needs, because a
|
||||
// user acts on exactly one geometry: the one the result was produced with.
|
||||
@@ -593,6 +680,64 @@ ReportDocument BuildReportDocument(const std::string &output_prefix,
|
||||
doc.sections.push_back(std::move(s));
|
||||
}
|
||||
|
||||
// ------------------------------------------------- XDS GEOMETRY KEYWORDS
|
||||
// The geometry the result was integrated at - on a rotation run the post-refined one - as the XDS.INP
|
||||
// keywords that carry it, so it can be pasted into an XDS run rather than retyped. Geometry only:
|
||||
// XDS determines the cell and the space group itself. Comment lines, like all prose here, with a
|
||||
// fixed prefix so `sed -n 's/^# XDS.INP> //p'` lifts them.
|
||||
//
|
||||
// Following the XDS conventions, Kabsch (2010) Acta Cryst. D66, 125-132: XDS places a pixel at
|
||||
// x_lab(i,j) = (i-ORGX)*QX*X_axis + (j-ORGY)*QY*Y_axis + DETECTOR_DISTANCE*(X_axis x Y_axis)
|
||||
// in a laboratory frame that is ours (x along a column step, y along a row step, z along the beam),
|
||||
// with pixels numbered from 1 and ORGX/ORGY the foot of the perpendicular from the crystal - our PONI
|
||||
// (docs/DETECTOR_GEOMETRY.md). A tilt therefore goes out as the two detector axes, taken as
|
||||
// differences of LabCoord.
|
||||
//
|
||||
// ROTATION_AXIS is rugnux's NEGATED: rugnux turns the crystal by -phi about its axis (IndexAndRefine
|
||||
// takes the lattice to an image with GetTransformationAngle(-angle)), XDS by a right-handed +phi about
|
||||
// ROTATION_AXIS. Checked with IDXREF on two sets given this geometry: the negated axis indexes 99.5%
|
||||
// and 14% of the spots, rugnux's own sign 37% and 0.04% (the chance floor of a mis-signed spindle).
|
||||
// The axis is the one the run integrated at, after any sign rescue.
|
||||
if (std::isfinite(result.used_beam_x_pxl) && std::isfinite(result.used_distance_mm)) {
|
||||
constexpr double RAD = PI / 180.0;
|
||||
const auto &t = result.used_detector_tilt_deg;
|
||||
DiffractionExperiment x = experiment;
|
||||
x.BeamX_pxl(result.used_beam_x_pxl).BeamY_pxl(result.used_beam_y_pxl)
|
||||
.DetectorDistance_mm(result.used_distance_mm)
|
||||
.PoniRot1_rad(static_cast<float>(t[0] * RAD)).PoniRot2_rad(static_cast<float>(t[1] * RAD))
|
||||
.PoniRot3_rad(static_cast<float>(t[2] * RAD));
|
||||
const auto geom = x.GetDiffractionGeometry();
|
||||
const float bx = geom.GetBeamX_pxl(), by = geom.GetBeamY_pxl(), pixel_mm = geom.GetPixelSize_mm();
|
||||
const Coord origin = geom.LabCoord(bx, by);
|
||||
const Coord x_axis = ((geom.LabCoord(bx + 1.0f, by) - origin) / pixel_mm).Normalize();
|
||||
const Coord y_axis = ((geom.LabCoord(bx, by + 1.0f) - origin) / pixel_mm).Normalize();
|
||||
|
||||
std::string k;
|
||||
const auto line = [&k](const std::string &l) { k += " XDS.INP> " + l + "\n"; };
|
||||
line(fmt::format("ORGX= {:.2f} ORGY= {:.2f}", bx + 1.0f, by + 1.0f));
|
||||
line(fmt::format("DETECTOR_DISTANCE= {:.3f}", geom.GetDetectorDistance_mm()));
|
||||
if (t[0] != 0.0 || t[1] != 0.0 || t[2] != 0.0) {
|
||||
line(fmt::format("DIRECTION_OF_DETECTOR_X-AXIS= {:.6f} {:.6f} {:.6f}", x_axis.x, x_axis.y, x_axis.z));
|
||||
line(fmt::format("DIRECTION_OF_DETECTOR_Y-AXIS= {:.6f} {:.6f} {:.6f}", y_axis.x, y_axis.y, y_axis.z));
|
||||
}
|
||||
line(fmt::format("INCIDENT_BEAM_DIRECTION= 0 0 1 X-RAY_WAVELENGTH= {:.6f}", geom.GetWavelength_A()));
|
||||
if (const auto gonio = experiment.GetGoniometer()) {
|
||||
const Coord axis = -result.used_rotation_axis.value_or(gonio->GetAxis());
|
||||
// + 0.0f: no "-0.000000" from a negated zero.
|
||||
line(fmt::format("ROTATION_AXIS= {:.6f} {:.6f} {:.6f}", axis.x + 0.0f, axis.y + 0.0f, axis.z + 0.0f));
|
||||
}
|
||||
k.pop_back();
|
||||
|
||||
ReportSection s;
|
||||
s.title = "XDS GEOMETRY KEYWORDS";
|
||||
Add(s, Prose(" The geometry this result was integrated at, as XDS.INP keywords - ORGX/ORGY numbered from\n"
|
||||
" 1, ROTATION_AXIS in XDS's sense of rotation. Geometry only; XDS determines the cell and\n"
|
||||
" space group itself. Lift them with: sed -n 's/^# XDS.INP> //p' <this file>"));
|
||||
Add(s, Blank());
|
||||
Add(s, Prose(k));
|
||||
doc.sections.push_back(std::move(s));
|
||||
}
|
||||
|
||||
// ------------------------------------------------------------- 2. CRYSTAL
|
||||
{
|
||||
ReportSection s;
|
||||
@@ -1452,6 +1597,53 @@ ReportDocument BuildReportDocument(const std::string &output_prefix,
|
||||
Add(s, Blank());
|
||||
}
|
||||
|
||||
// ---- overloads and overlapping spots: two ways the detector or the pattern took reflections
|
||||
// the crystal gave, each with what to change at the next collection
|
||||
if (result.integrated_reflections > 0) {
|
||||
const double pct = 100.0 * result.overloaded_reflections / result.integrated_reflections;
|
||||
const double pct_low = result.integrated_reflections_low_res > 0
|
||||
? 100.0 * result.overloaded_reflections_low_res / result.integrated_reflections_low_res : 0.0;
|
||||
Add(s, KeyInt("OVERLOADED_REFLECTIONS", result.overloaded_reflections));
|
||||
Add(s, KeyReal("OVERLOADED_REFLECTIONS_PCT", pct, "{:.3f}"));
|
||||
Add(s, KeyReal("OVERLOADED_REFLECTIONS_LOW_RES_PCT", pct_low, "{:.3f}"));
|
||||
overload_row = fmt::format("{:.3f}% of the integrated reflections, {:.3f}% of those coarser than {:.0f} A",
|
||||
pct, pct_low, OVERLOAD_LOW_RES_A);
|
||||
if (pct >= OVERLOAD_WARNING_PCT || pct_low >= OVERLOAD_LOW_RES_WARNING_PCT)
|
||||
Warn(doc, PathologyCode::OVERLOADS, fmt::format(
|
||||
"Overloads: {} integrated reflections ({:.3f}%; {:.2f}% of those coarser than {:.0f} A) held "
|
||||
"a saturated pixel and are not measurements - the merge leaves them out, and they are the "
|
||||
"strongest reflections the crystal gave. Attenuate the beam or shorten the exposure (or "
|
||||
"collect a second, attenuated low-resolution sweep)",
|
||||
result.overloaded_reflections, pct, pct_low, OVERLOAD_LOW_RES_A));
|
||||
}
|
||||
if (result.predicted_reflections > 0
|
||||
&& experiment.GetBraggIntegrationSettings().GetOverlap() != OverlapMode::Off) {
|
||||
const double pct = 100.0 * static_cast<double>(result.overlapped_reflections)
|
||||
/ static_cast<double>(result.predicted_reflections);
|
||||
const double pct_major = 100.0 * static_cast<double>(result.overlapped_reflections_major)
|
||||
/ static_cast<double>(result.predicted_reflections);
|
||||
const double minpk = 100.0 * experiment.GetBraggIntegrationSettings().GetOverlapMinPeak();
|
||||
Add(s, KeyReal("OVERLAPPED_REFLECTIONS_PCT", pct, "{:.2f}"));
|
||||
Add(s, KeyReal("OVERLAPPED_REFLECTIONS_MAJOR_PCT", pct_major, "{:.2f}"));
|
||||
overlap_row = fmt::format("{:.2f}% of the predicted reflections share signal pixels, {:.2f}% keep under "
|
||||
"{:.0f}% of their profile", pct, pct_major, minpk);
|
||||
if (pct_major >= OVERLAP_WARNING_PCT)
|
||||
Warn(doc, PathologyCode::OVERLAPS, fmt::format(
|
||||
"Overlapping spots: {:.2f}% of the predicted reflections keep less than {:.0f}% of their "
|
||||
"profile as their own ({:.2f}% share signal pixels with a neighbour) - the shared pixels "
|
||||
"are left out of the fit, which costs precision. Move the detector back, slice the "
|
||||
"rotation finer or use a smaller beam", pct_major, minpk, pct));
|
||||
}
|
||||
if (result.integrated_reflections > 0 || result.predicted_reflections > 0) {
|
||||
Add(s, Blank());
|
||||
Add(s, Prose(" OVERLOADED_REFLECTIONS counts the integrated reflections (every frame's partial) with a\n"
|
||||
" saturated pixel in the signal disk, OVERLOADED_REFLECTIONS_LOW_RES_PCT their share among\n"
|
||||
" those coarser than 4 A, where they concentrate. OVERLAPPED_REFLECTIONS_PCT is the share\n"
|
||||
" of the predicted reflections whose signal disk a nearer neighbour takes pixels from,\n"
|
||||
" OVERLAPPED_REFLECTIONS_MAJOR_PCT the share left with less than --overlap-minpk of it."));
|
||||
Add(s, Blank());
|
||||
}
|
||||
|
||||
// ---- translational pseudo-symmetry
|
||||
// Printed immediately before the twinning block, and in the same section, because the two
|
||||
// interact: a pseudo-translation that is not a half-integer vector biases <|L|> upwards and
|
||||
@@ -2377,6 +2569,12 @@ ReportDocument BuildReportDocument(const std::string &output_prefix,
|
||||
row("Indexing choice", fmt::format("alternative indexing {} - one sweep is indexed consistently, "
|
||||
"so it matters only against other data: index against those "
|
||||
"with -z or --model", result.indexing_ambiguity_operators));
|
||||
if (!beam_centre_row.empty())
|
||||
checked_row("Beam centre", {PathologyCode::BEAM_CENTRE_OFF}, beam_centre_row);
|
||||
if (!overload_row.empty())
|
||||
checked_row("Overloads", {PathologyCode::OVERLOADS}, overload_row);
|
||||
if (!overlap_row.empty())
|
||||
checked_row("Overlaps", {PathologyCode::OVERLAPS}, overlap_row);
|
||||
if (!result.powder.rings_q_recipA.empty() || result.consensus_cell) {
|
||||
const auto &pw = result.powder;
|
||||
checked_row("Powder", {PathologyCode::POWDER_RINGS}, pw.rings_q_recipA.empty() ? std::string("no rings detected")
|
||||
|
||||
+11
-68
@@ -309,73 +309,6 @@ Rugnux::Rugnux(JFJochReader &reader, DiffractionExperiment experiment,
|
||||
}
|
||||
}
|
||||
|
||||
// Report the detector geometry the way XDS states it, so a run can be compared line by line against
|
||||
// the IDXREF.LP / CORRECT.LP of XDS on the same data. XDS is never given this geometry - the durin
|
||||
// plugin hands it image data only (plugin_get_header returns dimensions, bytes per pixel, pixel size
|
||||
// and frame count, nothing else) and XDS refines its own from XDS.INP - which is precisely why
|
||||
// having ours in the same convention is worth printing.
|
||||
//
|
||||
// The two laboratory frames coincide: x along increasing detector column, y along increasing row,
|
||||
// z along the beam. XDS places a pixel at
|
||||
// x_lab(i,j) = (i-ORGX)*QX*X_axis + (j-ORGY)*QY*Y_axis + DETECTOR_DISTANCE*(X_axis x Y_axis)
|
||||
// and DiffractionGeometry::LabCoord is that same expression with X_axis = poni_rot*(1,0,0) and
|
||||
// Y_axis = poni_rot*(0,1,0). So the tilt goes out as the two axis vectors rather than as angles,
|
||||
// which is also how XDS reports it after refinement. Taking the axes as differences of LabCoord
|
||||
// keeps them in step with whatever the geometry currently is, tilt included.
|
||||
//
|
||||
// ORGX/ORGY are 1-based - XDS counts pixels from 1 - and are the foot of the perpendicular from the
|
||||
// crystal rather than the direct beam, which is what our beam centre is too
|
||||
// (docs/DETECTOR_GEOMETRY.md), so they map across without correction.
|
||||
//
|
||||
// ROTATION_AXIS is printed once per pass, and it is not always the vector stored in the file: the
|
||||
// sign rescue further down adopts the opposite sign when the file's own indexes almost nothing, and
|
||||
// the following pass then prints what it adopted. On the corpus the first and last line of a run
|
||||
// disagree in sign on 12 of 115 rotation runs, so the LAST one - the axis the run actually
|
||||
// integrated at - is the one to read.
|
||||
//
|
||||
// That sign is not reliably the one XDS wants - but it is not a constant either, so do NOT negate
|
||||
// this line on principle. It is paired with the detector axis directions printed just above it and
|
||||
// has to be read together with them: two rotation datasets whose detector axes are both +x/+y wanted
|
||||
// the negation, while a third whose axes are negative wanted the sign as printed. The measure is the
|
||||
// fraction of COLSPOT spots IDXREF can index - 8.9% as printed against 92.6% negated, and 1060/10080
|
||||
// against 9624/10080; the other way round, 43.9% as printed against 13.5%.
|
||||
//
|
||||
// Read the HIGH member of each pair. That statistic has a chance floor which reaches about 15% on a
|
||||
// dense reciprocal lattice - a null with the frame numbers permuted scored 14.9% against a real
|
||||
// 15.4% - so a low value means "IDXREF gave up", not a measured near-zero: both wrong-sign runs
|
||||
// above stopped on INSUFFICIENT PERCENTAGE (< 50%) OF INDEXED REFLECTIONS, one of them after
|
||||
// refining the right cell and rejecting 9020 spots as too far from their ideal positions, which is
|
||||
// the signature of a mis-signed spindle. So try both signs and keep the one that indexes; a sign
|
||||
// carried over from another file is worth nothing.
|
||||
//
|
||||
// One of those two runs read a file our own converter had written, whose axis default was the
|
||||
// negative sign. It therefore says what XDS does with a given sign, which is what this comment is
|
||||
// about, and nothing about what a foreign file stores.
|
||||
static void LogXDSGeometry(const DiffractionExperiment &experiment) {
|
||||
Logger logger("Rugnux");
|
||||
|
||||
const auto geom = experiment.GetDiffractionGeometry();
|
||||
const float beam_x = geom.GetBeamX_pxl();
|
||||
const float beam_y = geom.GetBeamY_pxl();
|
||||
const float pixel_mm = geom.GetPixelSize_mm();
|
||||
|
||||
const Coord origin = geom.LabCoord(beam_x, beam_y);
|
||||
const Coord x_axis = ((geom.LabCoord(beam_x + 1.0f, beam_y) - origin) / pixel_mm).Normalize();
|
||||
const Coord y_axis = ((geom.LabCoord(beam_x, beam_y + 1.0f) - origin) / pixel_mm).Normalize();
|
||||
|
||||
logger.Info("XDS convention: ORGX= {:.2f} ORGY= {:.2f} DETECTOR_DISTANCE= {:.4f}",
|
||||
beam_x + 1.0f, beam_y + 1.0f, geom.GetDetectorDistance_mm());
|
||||
logger.Info("XDS convention: DIRECTION_OF_DETECTOR_X-AXIS= {:.6f} {:.6f} {:.6f}",
|
||||
x_axis.x, x_axis.y, x_axis.z);
|
||||
logger.Info("XDS convention: DIRECTION_OF_DETECTOR_Y-AXIS= {:.6f} {:.6f} {:.6f}",
|
||||
y_axis.x, y_axis.y, y_axis.z);
|
||||
logger.Info("XDS convention: INCIDENT_BEAM_DIRECTION= 0 0 1 X-RAY_WAVELENGTH= {:.6f} QX= QY= {:.6f}",
|
||||
geom.GetWavelength_A(), pixel_mm);
|
||||
if (const auto goniometer = experiment.GetGoniometer())
|
||||
logger.Info("XDS convention: ROTATION_AXIS= {:.6f} {:.6f} {:.6f}",
|
||||
goniometer->GetAxis().x, goniometer->GetAxis().y, goniometer->GetAxis().z);
|
||||
}
|
||||
|
||||
namespace rugnux_internal {
|
||||
// The FFT recovers a basis vector only up to fft_max_unit_cell_A: FFTIndexer sizes its projected
|
||||
// histogram from that value and the transform's last usable bin IS that length. So a reference cell
|
||||
@@ -434,8 +367,17 @@ ProcessResult Rugnux::Run(RugnuxObserver *observer) {
|
||||
// on the rotation two-pass that is under a third of what the run actually took, the pre-scan and the
|
||||
// first pass being invisible in it. Time the whole thing here, where every pass is inside.
|
||||
const auto run_start = std::chrono::steady_clock::now();
|
||||
DiffractionExperiment header_experiment = experiment_;
|
||||
ProcessResult result = RunAllPasses(observer);
|
||||
result.total_time_s = std::chrono::duration<double>(std::chrono::steady_clock::now() - run_start).count();
|
||||
// For the report's beam-centre warning: what the background measured, and what this crystal asks
|
||||
// of the centre - at the file's geometry, as the pre-scan's check asks it, but on the cell the run
|
||||
// determined rather than the long axis that check has to assume.
|
||||
JoinBeamCenterCapture();
|
||||
result.background_beam_centre = measured_beam_center_;
|
||||
if (result.consensus_cell)
|
||||
header_experiment.SetUnitCell(*result.consensus_cell);
|
||||
result.beam_centre_need_pxl = BeamCenterNeed_pxl(header_experiment);
|
||||
return result;
|
||||
}
|
||||
|
||||
@@ -778,7 +720,6 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b
|
||||
: calibration ? "powder calibration" : "azimuthal integration";
|
||||
logger.Info("Processing {} images (range {}-{}, stride {}) using {} threads [{}]",
|
||||
images_to_process, start_image, end_image, config_.stride, config_.nthreads, mode_name);
|
||||
LogXDSGeometry(experiment_);
|
||||
if (observer)
|
||||
observer->OnPhase(full ? "Full analysis" : calibration ? "Powder calibration" : "Azimuthal integration");
|
||||
|
||||
@@ -1066,6 +1007,8 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b
|
||||
result.used_direct_beam_x_pxl = direct.first;
|
||||
result.used_direct_beam_y_pxl = direct.second;
|
||||
}
|
||||
if (const auto g = experiment_.GetGoniometer())
|
||||
result.used_rotation_axis = g->GetAxis();
|
||||
if (census_future.valid()) {
|
||||
TimingMark mark(logger, "leftover-lattice census wait");
|
||||
result.leftover_lattices = census_future.get();
|
||||
|
||||
@@ -301,6 +301,9 @@ struct ProcessResult {
|
||||
// rotations, which is a geometry no pass ever ran at.
|
||||
std::array<double, 3> used_detector_tilt_deg{};
|
||||
double used_direct_beam_x_pxl = 0.0, used_direct_beam_y_pxl = 0.0;
|
||||
// The rotation axis the same pass integrated at - the sign rescue can adopt the opposite of the
|
||||
// file's. The report's XDS keyword block reads it. Empty on stills.
|
||||
std::optional<Coord> used_rotation_axis;
|
||||
std::optional<UnitCell> consensus_cell;
|
||||
// The centring `consensus_cell` is written in; empty when nothing named one (the stills path
|
||||
// averages per-image cells and has no lattice search behind them). A conventional cell says
|
||||
@@ -530,6 +533,24 @@ struct ProcessResult {
|
||||
// did not post-refine (stills, single pass, or the refinement did not converge).
|
||||
std::optional<PostRefineResult> post_refine;
|
||||
|
||||
// The beam centre the scattered background puts the direct beam at (the pre-scan's check), and
|
||||
// how far the centre may be wrong on this crystal (BeamCenterNeed_pxl, with the run's cell where
|
||||
// it has one). Report only: the file's centre is read off the experiment the report is given, the
|
||||
// refined one off used_direct_beam_*. Empty / NaN where the check did not run.
|
||||
std::optional<BeamCenterEstimate> background_beam_centre;
|
||||
float beam_centre_need_pxl = NAN;
|
||||
|
||||
// Overloads: integrated reflections (partial observations, every frame) a saturated pixel of the
|
||||
// signal disk made unmeasurable (Reflection::overloaded), over the whole range and over the
|
||||
// reflections coarser than OVERLOAD_LOW_RES_A, where they concentrate. Counted on the integration
|
||||
// outcome before the merge; zero where nothing was integrated.
|
||||
int64_t integrated_reflections = 0, overloaded_reflections = 0;
|
||||
int64_t integrated_reflections_low_res = 0, overloaded_reflections_low_res = 0;
|
||||
// Overlapping spots, as the integrator's overlap treatment measured them (BraggIntegrationCounts):
|
||||
// predicted reflections whose signal disk a nearer neighbour takes pixels from, and of those the
|
||||
// ones left with less than --overlap-minpk of it. Zero with --overlap off.
|
||||
uint64_t predicted_reflections = 0, overlapped_reflections = 0, overlapped_reflections_major = 0;
|
||||
|
||||
// Which pass this result is (1-based) and how many ran, plus one line saying which pass became the
|
||||
// canonical output and on what evidence. pass_count == 1 outside the rotation two-pass.
|
||||
int pass_number = 1;
|
||||
@@ -561,6 +582,10 @@ struct ProcessResult {
|
||||
std::optional<double> spindle_lost_unique_fraction;
|
||||
};
|
||||
|
||||
// Overloads are counted separately over the reflections coarser than this (A) - the Wilson-plot
|
||||
// convention's 4 A (CalcGlobalWilsonB), where the strong low-resolution reflections that saturate sit.
|
||||
constexpr float OVERLOAD_LOW_RES_A = 4.0f;
|
||||
|
||||
// The rotation two-pass quality guard: whether the refined-geometry pass (pass2) is worse than the
|
||||
// header-geometry pass (pass1), judged on their search merges. Returns why, in the words the run
|
||||
// reports it in, or an empty string when the refined pass stands.
|
||||
|
||||
@@ -127,7 +127,7 @@ void Rugnux::ImagePass(PipelineLocals &p) {
|
||||
// finishes (BraggIntegrationCounts). Reading an engine's counts synchronises its stream, so it is
|
||||
// done once per worker at the end rather than per image.
|
||||
std::atomic<uint64_t> bragg_predicted = 0, bragg_starved = 0, bragg_starved_nb = 0,
|
||||
bragg_fallback = 0;
|
||||
bragg_fallback = 0, bragg_overlapped = 0, bragg_overlapped_major = 0;
|
||||
|
||||
// Calibration by spots: the pooled spot list of the whole run. Ring clustering is O(n^2) in it and
|
||||
// the Hough circle centre O(n^3) in its first few hundred, so each image contributes a fair share of
|
||||
@@ -300,6 +300,8 @@ void Rugnux::ImagePass(PipelineLocals &p) {
|
||||
bragg_starved += bragg.bkg_starved;
|
||||
bragg_starved_nb += bragg.bkg_starved_by_neighbour;
|
||||
bragg_fallback += bragg.profile_fallback;
|
||||
bragg_overlapped += bragg.overlapped;
|
||||
bragg_overlapped_major += bragg.overlapped_major;
|
||||
};
|
||||
|
||||
if (observer)
|
||||
@@ -475,6 +477,16 @@ void Rugnux::ImagePass(PipelineLocals &p) {
|
||||
experiment_.GetBraggIntegrationSettings().GetR3(),
|
||||
bragg_predicted.load(), frac(bragg_starved), frac(bragg_starved_nb),
|
||||
frac(bragg_fallback));
|
||||
// Overlapping spots: signal disks a nearer neighbour's prediction takes pixels from, as the
|
||||
// overlap treatment decided ownership. For the report (OVERLAPPED_REFLECTIONS).
|
||||
result.predicted_reflections = bragg_predicted;
|
||||
result.overlapped_reflections = bragg_overlapped;
|
||||
result.overlapped_reflections_major = bragg_overlapped_major;
|
||||
if (experiment_.GetBraggIntegrationSettings().GetOverlap() != OverlapMode::Off)
|
||||
logger.Info("Overlapping spots: {:.2f}% of the predicted reflections share signal pixels with a "
|
||||
"nearer neighbour, {:.2f}% keep less than {:.0f}% of their profile as their own",
|
||||
frac(bragg_overlapped), frac(bragg_overlapped_major),
|
||||
100.0 * experiment_.GetBraggIntegrationSettings().GetOverlapMinPeak());
|
||||
}
|
||||
|
||||
// Every image failing is a total failure, not a run that produced nothing: it used to be reported
|
||||
|
||||
@@ -298,6 +298,32 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) {
|
||||
phase("Scaling and merging");
|
||||
Step(geometry_prepass ? "Scaling and post-refinement" : "Scaling and merging");
|
||||
|
||||
// Overloads, for the report: the integrated reflections a saturated signal pixel made unmeasurable.
|
||||
// Counted before anything below drops or flags one. The merge leaves them out (a rotation merge
|
||||
// drops the whole rocking event) - they are the strongest reflections, mostly at low resolution,
|
||||
// and what a model's low-resolution terms rest on.
|
||||
{
|
||||
for (const auto &io : indexer->GetIntegrationOutcome())
|
||||
for (const auto &r : io.reflections) {
|
||||
const bool low_res = r.d >= OVERLOAD_LOW_RES_A;
|
||||
++result.integrated_reflections;
|
||||
result.integrated_reflections_low_res += low_res;
|
||||
if (r.overloaded) {
|
||||
++result.overloaded_reflections;
|
||||
result.overloaded_reflections_low_res += low_res;
|
||||
}
|
||||
}
|
||||
if (result.integrated_reflections > 0)
|
||||
logger.Info("Overloads: {} of {} integrated reflections ({:.3f}%) held a saturated pixel, {:.3f}% of "
|
||||
"the {} coarser than {:.0f} A", result.overloaded_reflections,
|
||||
result.integrated_reflections,
|
||||
100.0 * result.overloaded_reflections / result.integrated_reflections,
|
||||
result.integrated_reflections_low_res > 0
|
||||
? 100.0 * result.overloaded_reflections_low_res / result.integrated_reflections_low_res
|
||||
: 0.0,
|
||||
result.integrated_reflections_low_res, OVERLOAD_LOW_RES_A);
|
||||
}
|
||||
|
||||
// Ice-ring handling (--detect-ice-rings): flag reflections sitting on a hexagonal-ice powder
|
||||
// ring. Their integrated intensity is contaminated by the strong, variable ice background, so
|
||||
// they drag the per-image scale fit. Flag them here so scaling (the per-image G fit and the
|
||||
|
||||
@@ -1200,3 +1200,96 @@ TEST_CASE("ResultReport_Checks", "[Diagnostics]") {
|
||||
CHECK(std::ranges::any_of(c.codes, [&](const std::string &code) {
|
||||
return std::ranges::find(doc.pathology_flags, code) != doc.pathology_flags.end(); }));
|
||||
}
|
||||
|
||||
TEST_CASE("ResultReport_BeamCentreOff", "[Diagnostics]") {
|
||||
DiffractionExperiment x(DetJF(1));
|
||||
x.BeamX_pxl(500.0f).BeamY_pxl(400.0f);
|
||||
const auto file = x.GetDiffractionGeometry().GetDirectBeam_pxl();
|
||||
ProcessResult result;
|
||||
result.post_refine = PostRefineResult{};
|
||||
result.beam_centre_need_pxl = 2.0f;
|
||||
result.background_beam_centre = BeamCenterEstimate{.beam_x_pxl = static_cast<float>(file.first) + 7.0f,
|
||||
.beam_y_pxl = static_cast<float>(file.second),
|
||||
.sigma_pxl = 0.5f};
|
||||
|
||||
// Refined within the tolerance: the background's word is not enough where the centre was refined.
|
||||
result.used_direct_beam_x_pxl = file.first + 1.0;
|
||||
result.used_direct_beam_y_pxl = file.second;
|
||||
auto text = RenderResultReport("p", "in.h5", x, result);
|
||||
CHECK(text.find("BEAM_CENTRE_OFF") == std::string::npos);
|
||||
CHECK(text.find("\nBEAM_CENTRE_SHIFT= 1.00\n") != std::string::npos);
|
||||
CHECK(text.find("\nBEAM_CENTRE_TOLERANCE= 2.00\n") != std::string::npos);
|
||||
|
||||
// 4 px off: over the tolerance but under three times it - stated, not warned.
|
||||
result.used_direct_beam_x_pxl = file.first;
|
||||
result.used_direct_beam_y_pxl = file.second + 4.0;
|
||||
text = RenderResultReport("p", "in.h5", x, result);
|
||||
CHECK(text.find("BEAM_CENTRE_OFF") == std::string::npos);
|
||||
CHECK(text.find("under three times it") != std::string::npos);
|
||||
|
||||
// Refined 8 px off the file's centre: warned, with both centres, px and mm, and the advice.
|
||||
result.used_direct_beam_y_pxl = file.second + 8.0;
|
||||
text = RenderResultReport("p", "in.h5", x, result);
|
||||
CHECK(text.find("\nPATHOLOGY_FLAGS= BEAM_CENTRE_OFF\n") != std::string::npos);
|
||||
CHECK(text.find("8.0 px (0.60 mm)") != std::string::npos);
|
||||
CHECK(text.find("correct the beam centre at the beamline") != std::string::npos);
|
||||
|
||||
// Nothing refined (stills): the background alone warns, worded as a check.
|
||||
result.post_refine.reset();
|
||||
result.used_direct_beam_x_pxl = file.first;
|
||||
result.used_direct_beam_y_pxl = file.second;
|
||||
text = RenderResultReport("p", "in.h5", x, result);
|
||||
CHECK(text.find("\nPATHOLOGY_FLAGS= BEAM_CENTRE_OFF\n") != std::string::npos);
|
||||
CHECK(text.find("the background's word alone") != std::string::npos);
|
||||
}
|
||||
|
||||
TEST_CASE("ResultReport_OverloadsAndOverlaps", "[Diagnostics]") {
|
||||
DiffractionExperiment x(DetJF(1));
|
||||
ProcessResult result;
|
||||
result.integrated_reflections = 100000;
|
||||
result.integrated_reflections_low_res = 5000;
|
||||
result.overloaded_reflections = 5;
|
||||
result.overloaded_reflections_low_res = 5;
|
||||
result.predicted_reflections = 100000;
|
||||
result.overlapped_reflections = 5000;
|
||||
result.overlapped_reflections_major = 500;
|
||||
|
||||
auto text = RenderResultReport("p", "in.h5", x, result);
|
||||
CHECK(text.find("\nOVERLOADED_REFLECTIONS= 5\n") != std::string::npos);
|
||||
CHECK(text.find("\nOVERLOADED_REFLECTIONS_LOW_RES_PCT= 0.100\n") != std::string::npos);
|
||||
CHECK(text.find("\nOVERLAPPED_REFLECTIONS_MAJOR_PCT= 0.50\n") != std::string::npos);
|
||||
CHECK(text.find("\nPATHOLOGY_FLAGS= NONE\n") != std::string::npos);
|
||||
|
||||
// 1.2% of the low-resolution reflections overloaded, 1.5% left under MINPK: both warn.
|
||||
result.overloaded_reflections = 60;
|
||||
result.overloaded_reflections_low_res = 60;
|
||||
result.overlapped_reflections_major = 1500;
|
||||
text = RenderResultReport("p", "in.h5", x, result);
|
||||
CHECK(text.find("\nPATHOLOGY_FLAGS= OVERLOADS OVERLAPS\n") != std::string::npos);
|
||||
CHECK(text.find("Attenuate the beam") != std::string::npos);
|
||||
CHECK(text.find("Move the detector back") != std::string::npos);
|
||||
}
|
||||
|
||||
TEST_CASE("ResultReport_XDSGeometryKeywords", "[Diagnostics]") {
|
||||
DiffractionExperiment x(DetJF(1));
|
||||
x.BeamX_pxl(500.0f).BeamY_pxl(400.0f).DetectorDistance_mm(100.0f);
|
||||
x.Goniometer(GoniometerAxis("omega", 10.0f, 0.1f, Coord(1, 0, 0), {}));
|
||||
ProcessResult result;
|
||||
result.consensus_cell = UnitCell{.a = 79.0f, .b = 79.0f, .c = 38.0f,
|
||||
.alpha = 90.0f, .beta = 90.0f, .gamma = 90.0f};
|
||||
// The refined geometry, not the file's, goes out - 1-based, untilted here.
|
||||
result.used_beam_x_pxl = 510.0f;
|
||||
result.used_beam_y_pxl = 395.5f;
|
||||
result.used_distance_mm = 120.0f;
|
||||
result.used_rotation_axis = Coord(-1, 0, 0);
|
||||
|
||||
const auto text = RenderResultReport("p", "in.h5", x, result);
|
||||
CHECK(text.find("XDS GEOMETRY KEYWORDS") != std::string::npos);
|
||||
CHECK(text.find("\n# XDS.INP> ORGX= 511.00 ORGY= 396.50\n") != std::string::npos);
|
||||
CHECK(text.find("\n# XDS.INP> DETECTOR_DISTANCE= 120.000\n") != std::string::npos);
|
||||
// XDS turns the crystal the other way about the axis it is given.
|
||||
CHECK(text.find("\n# XDS.INP> ROTATION_AXIS= 1.000000 0.000000 0.000000\n") != std::string::npos);
|
||||
// Untilted: no detector axes; and never the cell or the space group.
|
||||
CHECK(text.find("XDS.INP> DIRECTION_OF_DETECTOR") == std::string::npos);
|
||||
CHECK(text.find("XDS.INP> UNIT_CELL") == std::string::npos);
|
||||
}
|
||||
@@ -21,6 +21,7 @@ Left figure: R-free of the deposited model after the battery's REFMAC check on t
|
||||
(x) and on Rugnux's data (y). By default the refined protocol (model_check.py: the same model
|
||||
refined against each data set at its own resolution, R-free over the free reflections both data sets
|
||||
share); --rfree first-cycle draws the first-cycle rigid-body numbers at the common resolution instead.
|
||||
A run without REFMAC values (the check is opt-in, --model-check) gets no R-free figure; the rest are drawn.
|
||||
Right figure: the deposited high-resolution limit (x) against Rugnux's own cut
|
||||
(y), log axes so the low-resolution entries show. Small molecules (no deposited model) and sets
|
||||
whose verdict is `fail` (a wrong lattice: nothing to compare) are left out.
|
||||
@@ -75,7 +76,10 @@ set style line 2 pt 7 ps 0.35 lc rgb "#0072b2"
|
||||
# counts: names left-aligned, numbers right-aligned in a column beside them
|
||||
set label 11 "Structures\nFacilities\nBeamlines" at graph 0.05, graph 0.93 left front
|
||||
set label 12 "%(n)d\n%(nf)d\n%(nb)d" at graph 0.31, graph 0.93 right front
|
||||
"""
|
||||
|
||||
# Drawn only where the runs carry REFMAC R-free values; the battery's REFMAC check is opt-in.
|
||||
GNUPLOT_RFREE = r"""
|
||||
# --- R-free
|
||||
set output "%(out)s/rfree.pdf"
|
||||
set xlabel "{/:Italic R}_{free}, depositor data"
|
||||
@@ -86,7 +90,9 @@ set format x "%%.1f"; set format y "%%.1f"
|
||||
set label 2 "Within ±%(band).2f: %(within)d/%(n)d" at graph 0.95, graph 0.07 right front
|
||||
plot x ls 1, "%(out)s/rfree.dat" using 1:2 ls 2
|
||||
unset label 2
|
||||
"""
|
||||
|
||||
GNUPLOT_REST = r"""
|
||||
# --- high-resolution limit
|
||||
set output "%(out)s/dmin.pdf"
|
||||
set logscale xy
|
||||
@@ -255,7 +261,8 @@ def main():
|
||||
else "set key off")
|
||||
params = dict(out=out, tkey=tkey, n=n, nf=nf, nb=nb, dn=dn, dnf=dnf, dnb=dnb, within=within, band=RFREE_BAND,
|
||||
deeper=deeper, shallower=sum(y > x for x, y, _ in dmin),
|
||||
rlo=0.05 * math.floor((min(values) - 0.01) / 0.05), rhi=0.05 * math.ceil((max(values) + 0.01) / 0.05),
|
||||
rlo=0.05 * math.floor((min(values, default=0.2) - 0.01) / 0.05),
|
||||
rhi=0.05 * math.ceil((max(values, default=0.3) + 0.01) / 0.05),
|
||||
dlo=0.5, dhi=8.0,
|
||||
bin=TIME_BIN, thi=TIME_BIN * math.ceil(max(secs) / TIME_BIN), tyhi=10 * math.ceil((tallest + 1) / 10),
|
||||
tmed=tmed, tn=len(times),
|
||||
@@ -263,9 +270,9 @@ def main():
|
||||
zyhi=5 * math.ceil((max(collections.Counter(math.floor(t / ZOOM_BIN) for t in secs
|
||||
if t < ZOOM_MAX).values()) + 1) / 5))
|
||||
gp = os.path.join(out, "figures.gp")
|
||||
open(gp, "w").write(GNUPLOT % params)
|
||||
open(gp, "w").write((GNUPLOT + (GNUPLOT_RFREE if rfree else "") + GNUPLOT_REST) % params)
|
||||
subprocess.run(["gnuplot", gp], check=True)
|
||||
for fig in ("rfree", "dmin", "time", "time_zoom"): # the PNG is the PDF rasterised at 600 d.p.i.
|
||||
for fig in (["rfree"] if rfree else []) + ["dmin", "time", "time_zoom"]: # the PNG is the PDF rasterised at 600 d.p.i.
|
||||
subprocess.run(["pdftoppm", "-r", "600", "-png", "-singlefile", os.path.join(out, fig + ".pdf"),
|
||||
os.path.join(out, fig)], check=True)
|
||||
|
||||
@@ -276,8 +283,10 @@ def main():
|
||||
fh.write(f"run {path}\n label {m.get('label')} binary git {b.get('source_head')} "
|
||||
f"dirty {b.get('source_dirty')} sha256 {b.get('sha256')}\n flags {b.get('flags')} "
|
||||
f"runner git {m.get('runner_git')} dirty {m.get('runner_dirty')}\n")
|
||||
fh.write(f"R-free: {n} structures, {nf} facilities, {nb} beamlines; within +-{RFREE_BAND}: {within}\n"
|
||||
f"d_min: {dn} structures, {dnf} facilities, {dnb} beamlines; Rugnux reaches further on "
|
||||
fh.write((f"R-free: {n} structures, {nf} facilities, {nb} beamlines; within +-{RFREE_BAND}: {within}\n"
|
||||
if rfree else "R-free: no REFMAC values in these runs (run the battery with --model-check); "
|
||||
"no rfree figure\n")
|
||||
+ f"d_min: {dn} structures, {dnf} facilities, {dnb} beamlines; Rugnux reaches further on "
|
||||
f"{deeper}\n")
|
||||
queue_wait = max(time_rows[s]["elapsed_s"] - t for s, _, t, _, _ in times)
|
||||
fh.write(f"time: {len(times)} sets of arm(s) {', '.join(a.time_arms)}; field rugnux_wall_s "
|
||||
|
||||
Reference in new issue
Block a user