v1.0.0-rc.167 (#77)

* `rugnux --model` reports CC(model, data) - the correlation of the merged intensities with the placed, scaled model - by resolution shell, on the same shells as CC1/2, with the reflection count and a significance for each.
* `rugnux --model` fits the model's scale, anisotropic B and bulk-solvent parameters on the working reflections only, so the R-free it reports is measured against a model no free reflection helped scale.
* The bulk-solvent parameters of `rugnux --model` are searched over their physically meaningful range instead of being fitted without bounds, so a model is never scaled with a solvent term that has silently switched itself off.
* The rigid-body placement of `rugnux --model` uses the same bounded bulk solvent as the reported fit, so a model is no longer placed against a target carrying a solvent term with no physical meaning.
* `rugnux --model` puts the model into the data's own description of the lattice before placing it, so a model whose cell is written on other axes - I-centred where the run indexed C-centred, a different unique axis, a permuted orthorhombic cell - is placed rather than scored where it was read; `MODEL_CHANGE_OF_BASIS=` and `MODEL_SETTING_AS_READ=` report it when it happens.
* The rugnux results report opens with a summary - `VERDICT=` (`OK`, `WARNINGS`, `UNUSABLE`, `FAILED`), `VERDICT_TEXT=`, `PATHOLOGY_FLAGS=` with one closed-vocabulary code per condition that warned, and the `WARNING:` lines, which used to close the file - and the sections after it are renumbered 1-5 with no gaps.
* `rugnux --developer` writes the full results report - the pipeline-internal keys and the long explanations the default report now leaves out - and `--finalist-ledger` adds the evidence for every space group the search considered, not only the one it adopted.
* The results report warns when the merged data carry no usable signal and when too little of reciprocal space was measured inside the fitted resolution, and omits `FITTED_RESOLUTION` where the CC1/2 curve it is fitted on never falls off.
* rugnux detects translational pseudo-symmetry and reports it under the `PSEUDO_TRANSLATION` flag as `TNCS_DETECTED=` and the `TNCS_*` keys - a translation the merged data are exactly invariant under is reported as `UNDECLARED_LATTICE_TRANSLATION=` under `LATTICE_TRANSLATION` instead - and a detected pseudo-translation can no longer buy a false screw axis in the space-group search or hide a twin from the L-test (`L_TEST_VS_TNCS=`).
* The space-group search determines glide planes from zonal systematic absences, so a non-Sohncke space group such as P 2_1/c or Pbca is named where the run previously stopped at its Sohncke subgroup; `SOHNCKE_SPACE_GROUP=` carries the best Sohncke group beside it on every run that searched, and a centre of symmetry is never claimed.
* Where the cell metric carries more rotational symmetry than the Bravais class the indexer named, the extra rotations are put to the intensities and the space-group search is asked again on the metric's own cell - adopted only where the intensities confirm the higher symmetry - so a lattice that is nearly but not exactly hexagonal, or whose reduction landed in a sub-cell, still reaches its true point group.
* Systematic-absence calls rest on the evidence rather than on counts: a screw axis whose absent class the data show extinct is no longer refused because a handful of reflections in it read as present, and `SPACE_GROUP_ALTERNATIVES=` no longer drops a candidate that differs only on a zone the sweep never measured.
* A reference correlation measured on too few reflections is refused instead of scored zero, so a run given a reference MTZ is no longer reindexed on an operator that mapped almost everything outside the reference's coverage.
* A frame counts as indexed from 6 spots on its lattice rather than 9, so a weakly diffracting crystal whose frames cannot carry 9 is no longer refused the lattice it fits; `--min-indexed-spots` overrides it.
* `-C` accepts a known cell in any equivalent description - conventional or primitive, centred or not - instead of only the reduced primitive form, so a centred cell given the way it is published no longer makes the run report that it found no lattice.
* Each reflection is corrected for the sensor's quantum efficiency at the angle it meets the detector (attenuation lengths from the NIST tables, which also fixes the spot-width parallax term on CdTe) and for the attenuation of the flight path between the sample and its pixel; `--flight-path air|helium|vacuum` declares the medium - default air, since no file states it - and the report says what was assumed and what it was worth. The unmerged MTZ records the factors in new `QE` and `FLIGHT` columns beside `LP`, so raw counts are `I / LP * QE * FLIGHT`, and `_process.h5` in new optional `qe` and `flight` datasets.
* Rotation geometry post-refinement fits the crystal and the detector at once, against the observed spot positions and the observed rocking angles together, so the refined distance depends far less on how wrong the file's distance was.
* A coarsely sliced sweep integrates correctly: partials are joined into one rocking event by angle rather than by frame count, so two crossings of the Ewald sphere are no longer summed into one full, and at 0.5 degrees per image or coarser the per-frame geometry refinement accepts a spot whose miss the exposure's own rotation accounts for.
* `rugnux --mode scale` reports the detector tilt and direct beam of the geometry it re-scaled at, instead of zeros that read as a flat detector, and no longer warns that no image was indexed on a run whose lattice came from its input file.
* Every rotation run that determined a space group and merged reports what the mounting cost: `SPINDLE_LOST_UNIQUE_FRACTION=` is the fraction (0-1) of unique reflections the mounting made unmeasurable under the measured point group, also written to the master as `/entry/MX/spindleLostUniqueFraction` and what the mounting warning fires on; `SPINDLE_SYMMETRY_AXIS_ANGLE_DEG=` / `SPINDLE_SYMMETRY_AXIS_ORDER=` describe the mounting in the `--developer` report.
* Stills and grid scans carry a per-image `spindle_blind_fraction` - how much of a rotation sweep's blind cone this orientation would make unrecoverable, 0.5 and above calling for a second orientation - through the CBOR stream, HDF5 (`/entry/MX/spindleBlindFraction`), the plot and scan-result APIs, and the viewer and frontend plots; an absent value means the frame could not be assessed and is not a 0.
* The results report's `REPORT_VERSION` is 7.

Reviewed-on: #77
Co-authored-by: Filip Leonarski <filip.leonarski@psi.ch>

image_analysis/WriteReflections.cpp

@@ -598,7 +598,7 @@ float DetectorY(const Reflection &r) { return std::isfinite(r.observed_y) ? r.ob
 
 // Sum each rocking event into one full observation. A rotation reflection is integrated image by
 // image, so it arrives here as a run of partials over consecutive frames; the run is cut where the
-// 3D combine cuts it - same raw hkl, frames no further apart than MAX_FRAME_GAP - so the exported
+// 3D combine cuts it - same raw hkl, frames no further apart than RockingEventFrameGap - so the exported
 // file and rugnux's own merge see exactly the same events.
 // The parts are added, plainly, with their variances in quadrature, which is what every other
 // rotation program writes as a full. Nothing is divided by the partiality: FRACTIONCALC carries the
@@ -630,8 +630,9 @@ float DetectorY(const Reflection &r) { return std::isfinite(r.observed_y) ? r.ob
 // mosaicity, RotationScaleMerge::SmoothMosaicity), so clamping it to 1 would hide the spread and buy
 // a reading program nothing.
 std::vector<Reflection> SumRockingEvents(const std::vector<IntegrationOutcome> &outcomes,
-                                         double min_partiality, double min_captured_fraction) {
-    constexpr float MAX_FRAME_GAP = 2.0f;   // == RotationScaleMerge's: what makes one rocking event
+                                         double min_partiality, double min_captured_fraction,
+                                         float wedge_deg) {
+    const float max_frame_gap = RockingEventFrameGap(wedge_deg);   // == RotationScaleMerge's
 
     // The sort key travels with the part instead of being read back through the pointer, the way the
     // merge's own ingest sort carries it (RotationScaleMerge's SortKey): there are millions of parts
@@ -661,20 +662,23 @@ std::vector<Reflection> SumRockingEvents(const std::vector<IntegrationOutcome> &
         size_t j = i + 1;
         while (j < parts.size() && parts[j].h == parts[i].h && parts[j].k == parts[i].k
                && parts[j].l == parts[i].l
-               && parts[j].image_number - parts[j - 1].image_number <= MAX_FRAME_GAP)
+               && parts[j].image_number - parts[j - 1].image_number <= max_frame_gap)
             ++j;
 
         double sum_p = 0.0, sum_I = 0.0, sum_var = 0.0, sum_var_bkg = 0.0;
-        double p_rlp = 0.0, p_frame = 0.0, p_x = 0.0, p_y = 0.0, p_delta_phi = 0.0, p_zeta = 0.0,
-               p_bkg = 0.0;
+        double p_corr = 0.0, p_qe = 0.0, p_flight = 0.0, p_frame = 0.0, p_x = 0.0, p_y = 0.0,
+               p_delta_phi = 0.0, p_zeta = 0.0, p_bkg = 0.0;
         for (size_t m = i; m < j; ++m) {
             const Reflection &r = *parts[m].r;
             const double p = r.partiality;
+            const float pc = r.prescaling_corr * r.qe_corr * r.flight_corr;   // LP x QE x flight path
             sum_p += p;
-            sum_I += static_cast<double>(r.I) * r.rlp;
-            sum_var += static_cast<double>(r.sigma) * r.sigma * r.rlp * r.rlp;
-            sum_var_bkg += static_cast<double>(r.var_bkg) * r.rlp * r.rlp;
-            p_rlp += p * r.rlp;
+            sum_I += static_cast<double>(r.I) * pc;
+            sum_var += static_cast<double>(r.sigma) * r.sigma * pc * pc;
+            sum_var_bkg += static_cast<double>(r.var_bkg) * pc * pc;
+            p_corr += p * pc;
+            p_qe += p * r.qe_corr;
+            p_flight += p * r.flight_corr;
             p_frame += p * r.image_number;
             p_x += p * DetectorX(r);
             p_y += p * DetectorY(r);
@@ -687,13 +691,18 @@ std::vector<Reflection> SumRockingEvents(const std::vector<IntegrationOutcome> &
         if (sum_p < min_partiality || sum_p < min_captured_fraction)
             continue;
 
-        // The Lorentz-polarization factor is applied by the writer, which multiplies I by rlp, so
-        // divide the event's own factor back out of the sums here. LP is the same geometry for every
-        // part of one event to a median 2e-4, so the file's I/LP is still the raw count sum.
-        full.rlp = static_cast<float>(p_rlp / sum_p);
-        full.I = static_cast<float>(sum_I / full.rlp);
-        full.sigma = static_cast<float>(std::sqrt(sum_var) / full.rlp);
-        full.var_bkg = static_cast<float>(sum_var_bkg / (static_cast<double>(full.rlp) * full.rlp));
+        // The prescaling factor is applied by the writer, which multiplies I by the whole correction, so
+        // divide the event's own factor back out of the sums here. It is the same geometry for every part
+        // of one event to a median 2e-4, so the file's I / LP * QE * FLIGHT is still the raw count sum. The
+        // event's mean is taken on the whole correction and on the two path terms; the
+        // Lorentz-polarization part is what is left when those are divided out of it.
+        const float mean_corr = static_cast<float>(p_corr / sum_p);
+        full.qe_corr = static_cast<float>(p_qe / sum_p);
+        full.flight_corr = static_cast<float>(p_flight / sum_p);
+        full.prescaling_corr = mean_corr / (full.qe_corr * full.flight_corr);
+        full.I = static_cast<float>(sum_I / mean_corr);
+        full.sigma = static_cast<float>(std::sqrt(sum_var) / mean_corr);
+        full.var_bkg = static_cast<float>(sum_var_bkg / (static_cast<double>(mean_corr) * mean_corr));
         full.partiality = static_cast<float>(sum_p);
         full.image_number = static_cast<float>(p_frame / sum_p);
         full.observed_x = static_cast<float>(p_x / sum_p);
@@ -724,9 +733,23 @@ void WriteUnmergedMtzReflections(const std::vector<IntegrationOutcome> &outcomes
     ds.wavelength = experiment.GetWavelength_A();
 
     // Every data column belongs to HKL_base and only the batches carry a dataset id, as in a
-    // POINTLESS file. Down to FLAG this is POINTLESS's own column set; the four after it are what
-    // rugnux measures beyond it - the offset of the reflection from the centre of its rocking curve,
-    // the Lorentz geometry of that curve, and the background that was subtracted.
+    // POINTLESS file. Down to FLAG this is POINTLESS's own column set, plus QE; the four after it
+    // are what rugnux measures beyond it - the offset of the reflection from the centre of its
+    // rocking curve, the Lorentz geometry of that curve, and the background that was subtracted.
+    //
+    // QE is not POINTLESS's, it is DIALS's: the detector's angle-dependent quantum efficiency, kept
+    // OUT of LP and given its own column because that is what both reference implementations do.
+    // LP here is Lorentz x polarization and nothing else, as XDS's RLP is (measured: XDS's column is
+    // flat to 0.1% across a detector over which the efficiency term spans 7%) and as DIALS's is
+    // (its LP comes from Corrections::lp(), lorentz/polarization, and it writes QE separately -
+    // a column of ones when it has no correction to report). Folding the sensor response into LP
+    // made this file assert something false about itself: a program that divided LP out to recover
+    // raw counts stripped the efficiency with it, and on an ordinary 13 keV geometry that term
+    // spans 13% end to end.
+    // Direction: QE is a DIVISOR, as in DIALS, and is normalised to 1 at normal incidence, so it
+    // rises above 1 toward the detector edge where a longer crossing makes the sensor more
+    // efficient. DIALS stores the un-normalised absorbed fraction instead; the two differ by the
+    // per-dataset constant QE(0), which is an overall scale.
     mtz.add_column("M/ISYM", 'Y', 0, -1, false);
     mtz.add_column("BATCH", 'B', 0, -1, false);
     mtz.add_column("I", 'J', 0, -1, false);
@@ -736,6 +759,8 @@ void WriteUnmergedMtzReflections(const std::vector<IntegrationOutcome> &outcomes
     mtz.add_column("YDET", 'R', 0, -1, false);
     mtz.add_column("ROT", 'R', 0, -1, false);
     mtz.add_column("LP", 'R', 0, -1, false);
+    mtz.add_column("QE", 'R', 0, -1, false);
+    mtz.add_column("FLIGHT", 'R', 0, -1, false);
     mtz.add_column("FLAG", 'I', 0, -1, false);
     mtz.add_column("DELPHI", 'R', 0, -1, false);
     mtz.add_column("ZETA", 'R', 0, -1, false);
@@ -761,11 +786,13 @@ void WriteUnmergedMtzReflections(const std::vector<IntegrationOutcome> &outcomes
     // them there, so the index the reflection was actually measured at is recoverable - that is the
     // crystal-frame information careless scales on, and what makes the file unmerged rather than a
     // merge waiting to happen.
-    // I and SIGI are the integrated intensity with the Lorentz-polarization factor applied and
-    // nothing else, which is what IOBS means in every unmerged format (LP records the factor, so the
-    // raw counts are I/LP). LP is geometry, not a scale, and a program that reads this file has no
-    // way to recover it. The partiality is NOT divided out - that is a scale, FRACTIONCALC carries
-    // it, and every program this file is for wants to handle it its own way.
+    // I and SIGI are the integrated intensity with every deterministic per-reflection correction
+    // applied and nothing else, which is what IOBS means in every unmerged format. Two columns
+    // record what was applied, so the raw count sum is recoverable exactly as I / LP * QE:
+    // LP is Lorentz x polarization and QE the sensor efficiency (a divisor, hence the multiply).
+    // Neither is a scale, and a program that reads this file has no way to recompute them. The
+    // partiality is NOT divided out - that IS a scale, FRACTIONCALC carries it, and every program
+    // this file is for wants to handle it its own way.
     gemmi::UnmergedHklMover hkl_mover(mtz.spacegroup);
     // The batch headers are written for every image the observations span, not only for the images
     // that produced one. AIMLESS starts a new run wherever the phi series jumps, so an image that
@@ -799,13 +826,20 @@ void WriteUnmergedMtzReflections(const std::vector<IntegrationOutcome> &outcomes
         mtz.data.push_back(static_cast<float>(hkl[2]));
         mtz.data.push_back(static_cast<float>((partials ? 256 : 0) + isym));
         mtz.data.push_back(static_cast<float>(batch));
-        mtz.data.push_back(r.I * r.rlp);
-        mtz.data.push_back(r.sigma * r.rlp);
+        const float corr = r.prescaling_corr * r.qe_corr * r.flight_corr;
+        mtz.data.push_back(r.I * corr);
+        mtz.data.push_back(r.sigma * corr);
         mtz.data.push_back(r.partiality);
         mtz.data.push_back(DetectorX(r));
         mtz.data.push_back(DetectorY(r));
         mtz.data.push_back(phi_start_deg(r.image_number) + wedge_deg / 2.0f);
-        mtz.data.push_back(r.rlp);
+        // prescaling_corr is Lorentz x polarization, which is what LP means; qe_corr is the sensor
+        // efficiency and flight_corr the flight path beside it, each written in the divide-by
+        // direction, so that I / LP * QE * FLIGHT is the raw count. FLIGHT is <= 1 where QE is >= 1:
+        // the sensor makes an oblique reflection read high and the medium makes it read low.
+        mtz.data.push_back(r.prescaling_corr);
+        mtz.data.push_back(1.0f / r.qe_corr);
+        mtz.data.push_back(1.0f / r.flight_corr);
         mtz.data.push_back(0.0f);            // FLAG: nothing here is a rejected observation
         mtz.data.push_back(r.delta_phi_deg);
         mtz.data.push_back(r.zeta);
@@ -815,7 +849,8 @@ void WriteUnmergedMtzReflections(const std::vector<IntegrationOutcome> &outcomes
     if (scanning && sum_partials) {
         for (const auto &r : SumRockingEvents(outcomes,
                                               experiment.GetScalingSettings().GetMinPartiality(),
-                                              experiment.GetScalingSettings().GetMinCapturedFraction()))
+                                              experiment.GetScalingSettings().GetMinCapturedFraction(),
+                                              wedge_deg))
             add_row(r);
     } else {
         for (const auto &outcome : outcomes)