rugnux: --search-min-zeta drops badly-measured observations from the symmetry search

zeta is the sine of the angle between a reflection's rocking path and the
spindle. Near 0 the reflection crosses the Ewald sphere almost tangentially,
spends many frames in diffracting position and is measured worst. The de-novo
space-group search asks how EQUAL an operator's paired intensities are, so its
answer is dominated by whichever reflections are measured worst - and when the
spindle lies in a lattice plane, an operator that permutes the two in-plane
axes samples a different mixture of measurement qualities than one that only
flips signs. That is not a fair comparison, and it can make a real symmetry
operator look like a twin law.

Measured on a thaumatin set mounted that way (its 4-fold is 88.9 deg from the
spindle), the added operators' disagreement is 1.74x the parent's over pairs
where both reflections have zeta < 0.85 and 1.003x - i.e. the symmetry is
exact - over pairs where both are above it. The search consequently refuses the
422 promotion and merges the crystal in P222, while the same data forced to the
right group give CC1/2 99.2% at multiplicity 10.7, matching XDS.

With the option the de-novo pass ignores those observations (the final merge
keeps everything - there completeness is the point):

  zeta cut   observations ignored   H ratio   adopted
  0 (off)                       -      1.47   P222
  0.5                     1620648      1.44   P222
  0.7                     3006013      1.34   P21212
  0.85                    4536724   promoted  P4212   (correct point group)

OFF BY DEFAULT, and it must stay off, because the same cut costs four other
crystals their space group (P41212 -> P212121, I23 -> P2, I23 -> I222 twice):
at 0.85 it discards 40-80% of all observations, which on a crystal whose
geometry is not the problem simply starves the search. Two independent
implementations - filtering the pairs that enter the statistic, and filtering
the observations that enter the merge - trade exactly the same crystals, so
this is a property of the cut and not of where it is applied. Verified
bit-identical to the previous binary when off.

The companion diagnostic is already there: the run now reports how close a
symmetry axis lies to the spindle, which is the geometry that makes this
option worth reaching for.

Implementation note for anyone tempted by the cheaper route: excluding these
observations from the ASU grouping alone does NOT work. The 3D combine selects
partials on corr, not on their group, so their intensity still reaches the
fulls and the merged intensities are unchanged - measured, the statistic did
not move by 0.03 while 67% of observations were nominally excluded. Zeroing
corr is what removes an observation from the combine, the merge and the error
model alike.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>

common/ScalingSettings.cpp

@@ -89,6 +89,18 @@ double ScalingSettings::GetMinCCForImage() const {
     return min_cc_for_image;
 }
 
+double ScalingSettings::GetSearchMinZeta() const {
+    return search_min_zeta;
+}
+
+ScalingSettings &ScalingSettings::SearchMinZeta(double input) {
+    if (input < 0.0 || input >= 1.0)
+        throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
+                              "Search zeta limit must be in [0,1)");
+    search_min_zeta = input;
+    return *this;
+}
+
 ScalingSettings &ScalingSettings::MinCCForImage(double input) {
     if (input < 0.0 || input > 1.0)
         throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "Min CC for image must be between 0 and 1");

common/ScalingSettings.h

@@ -28,6 +28,10 @@ class ScalingSettings {
     // min_partiality, which gates individual partials; this gates the assembled full. 0 = off (baseline).
     double min_captured_fraction = 0.0;
     double min_cc_for_image = 0.0;
+
+    // Exclude observations whose Lorentz geometry |zeta| falls below this from the DE-NOVO space-group
+    // search merge only (see RotationScaleMerge::search_min_zeta). 0 = off.
+    double search_min_zeta = 0.0;
     double outlier_reject_nsigma = 0.0; // per-observation merge outlier rejection (XDS/DIALS-style); 0 = off, e.g. 6 enables
 
     // Scale fulls: after the rotation 3D combine, refit a per-frame scale on the combined fulls (XDS
@@ -92,6 +96,7 @@ public:
     ScalingSettings& CaptureUncertaintyCoeff(double input);
     ScalingSettings& MinCapturedFraction(double input);
     ScalingSettings& MinCCForImage(double min_cc_for_image);
+    ScalingSettings& SearchMinZeta(double search_min_zeta);
     ScalingSettings& OutlierRejectNsigma(double input);
     ScalingSettings& ScaleFulls(bool input);
     ScalingSettings& AbsorptionIter(int input);
@@ -127,6 +132,7 @@ public:
     [[nodiscard]] double GetCaptureUncertaintyCoeff() const;
     [[nodiscard]] double GetMinCapturedFraction() const;
     [[nodiscard]] double GetMinCCForImage() const;
+    [[nodiscard]] double GetSearchMinZeta() const;
     [[nodiscard]] double GetOutlierRejectNsigma() const;
     [[nodiscard]] bool GetScaleFulls() const;
     [[nodiscard]] int GetAbsorptionIter() const;

image_analysis/scale_merge/RotationScaleMerge.cpp

@@ -177,6 +177,7 @@ RotationScaleMerge::RotationScaleMerge(const DiffractionExperiment &experiment,
     capture_uncertainty_coeff = s.GetCaptureUncertaintyCoeff();
     min_captured_fraction = s.GetMinCapturedFraction();
     min_cc_for_image = s.GetMinCCForImage();
+    search_min_zeta = s.GetSearchMinZeta();
     reject_nsigma = s.GetOutlierRejectNsigma();
     reject_outliers = reject_nsigma > 0.0;
     rfree_fraction = s.GetRfreeFraction();
@@ -1995,6 +1996,30 @@ RotationScaleMerge::Result RotationScaleMerge::Run(bool for_search,
     }
     FinalizePerFrameScale(cc, cc_n, partial_scaled);
 
+    // --- 2a. On the de-novo search pass only, drop observations whose Lorentz geometry is poor.
+    //         Zeroing corr is what removes an observation everywhere: the 3D combine, the merge and the
+    //         error model all require corr > 0. Excluding them from the ASU grouping alone is NOT
+    //         enough - the combine selects on corr, so their intensity would still reach the fulls. ---
+    if (for_search && search_min_zeta > 0.0) {
+        int64_t n_dropped = 0;
+        for (auto &o : partials)
+            if (!(std::isfinite(o.zeta) && o.zeta >= search_min_zeta)) {
+                if (std::isfinite(o.corr) && o.corr > 0.0f) ++n_dropped;
+                o.corr = 0.0f;
+            }
+#ifdef JFJOCH_USE_CUDA
+        if (gpu_active_ && n_dropped > 0) {
+            std::vector<float> corr(partials.size());
+            for (size_t i = 0; i < partials.size(); ++i) corr[i] = partials[i].corr;
+            gpu_->SetCorr(corr.data());
+        }
+#endif
+        if (n_dropped > 0)
+            logger.Info("Space-group search: ignoring {} observations with |zeta| < {:.2f} "
+                        "(they cross the Ewald sphere near-tangentially and are measured worst)",
+                        n_dropped, search_min_zeta);
+    }
+
     // --- 2b. Drop frames that do not agree with the merged reference (--min-image-cc). ---
     if (min_cc_for_image > 0.0) {
         std::vector<uint8_t> reject(n_frames, 0);

image_analysis/scale_merge/RotationScaleMerge.h

@@ -106,6 +106,16 @@ private:
     // indexes and still integrates - it just measures something that is not the crystal's diffraction,
     // and nothing downstream removes it. 0 = off.
     double min_cc_for_image = 0.0;
+
+    // Exclude observations with |zeta| below this from the DE-NOVO SEARCH merge only (the final merge
+    // keeps everything). zeta is the sine of the angle between a reflection's rocking path and the
+    // spindle: near 0 it crosses the Ewald sphere almost tangentially, spends many frames in
+    // diffracting position and is measured badly. The symmetry search compares how equal an operator's
+    // paired intensities are, so it is answered by whichever reflections are worst measured - and when
+    // the spindle lies in a lattice plane, an operator permuting the two in-plane axes samples a
+    // different mixture of qualities than one that only flips signs, which is not a fair comparison.
+    // Unlike a bound on I/sigma this is pure geometry, identical in meaning on every dataset. 0 = off.
+    double search_min_zeta = 0.0;
     double reject_nsigma = 0.0;
     bool reject_outliers = false;
     double rfree_fraction = 0.0;

rugnux/rugnux_cli.cpp

@@ -109,6 +109,7 @@ void print_usage() {
     std::cout << "   --min-captured-fraction <num>   rot3d: drop a combined full whose rocking curve was captured below this fraction (edge-of-sweep truncated fulls) (default: 0.7 for rotation, 0 otherwise; 0 = off)" << std::endl;
     std::cout << "   --mosaicity <num>               Diagnostic: fix the scaling mosaicity (deg) instead of the per-image seed" << std::endl;
     std::cout << "   --reject-outliers <num>         Per-observation merge outlier rejection, N sigma from the per-reflection median (default: 6 for rot3d, XDS/DIALS-style; 0 = off)" << std::endl;
+    std::cout << "   --search-min-zeta <num>         De-novo space-group search only: ignore observations whose Lorentz geometry |zeta| is below this (0-1). Reflections crossing the Ewald sphere near-tangentially are measured badly and can make a symmetry operator look wrong (default: 0 = use all)" << std::endl;
     std::cout << "   --min-image-cc <num>            Per-image CC limit in percent (default: no limit)" << std::endl;
     std::cout << "   --scaling-iterations <num>      Number of scaling iterations with no reference data (default: 3)" << std::endl;
     std::cout << "   -z, --reference-mtz <file>      Reference MTZ file" << std::endl;
@@ -152,6 +153,7 @@ enum {
     OPT_MAX_SPOTS,
     OPT_MIN_PARTIALITY,
     OPT_MIN_IMAGE_CC,
+    OPT_SEARCH_MIN_ZETA,
     OPT_SCALING_ITERATIONS,
     OPT_SCALING_HIGH_RESOLUTION,
     OPT_RESOLUTION_CUTOFF,
@@ -263,6 +265,7 @@ static option long_options[] = {
     {"min-captured-fraction", required_argument, nullptr, OPT_MIN_CAPTURED_FRACTION},
     {"mosaicity", required_argument, nullptr, OPT_MOSAICITY},
     {"min-image-cc", required_argument, nullptr, OPT_MIN_IMAGE_CC},
+    {"search-min-zeta", required_argument, nullptr, OPT_SEARCH_MIN_ZETA},
     {"scaling-iterations", required_argument, nullptr, OPT_SCALING_ITERATIONS},
     {"scaling-high-resolution", required_argument, nullptr, OPT_SCALING_HIGH_RESOLUTION},
     {"resolution-cutoff", required_argument, nullptr, OPT_RESOLUTION_CUTOFF},
@@ -527,6 +530,7 @@ int main(int argc, char **argv) {
     std::optional<double> capture_uncertainty_arg;   // explicit --capture-uncertainty; default depends on rot3d
     std::optional<double> forced_mosaicity_arg;      // diagnostic: fix the scaling mosaicity (deg) instead of the per-image seed
     double min_image_cc = 0.0;
+    double search_min_zeta = 0.0;   // --search-min-zeta; de-novo search merge only
     int64_t scaling_iter = 3;
     std::optional<CrystalLattice> forced_rotation_lattice;
     std::optional<int> refine_geometry;   // --refine-geometry[=N]: stills global geometry-refinement pass
@@ -851,6 +855,9 @@ int main(int argc, char **argv) {
         case OPT_MIN_IMAGE_CC:
             min_image_cc = parse_double_arg(optarg, "--min-image-cc", logger);
             break;
+        case OPT_SEARCH_MIN_ZETA:
+            search_min_zeta = parse_double_arg(optarg, "--search-min-zeta", logger);
+            break;
         case OPT_SCALING_HIGH_RESOLUTION:
             d_min_scale_merge = atof(optarg);
             break;
@@ -1425,6 +1432,7 @@ int main(int argc, char **argv) {
     scaling_settings.CaptureUncertaintyCoeff(capture_uncertainty_arg.value_or(rotation_indexing ? 1.0 : 0.0));
     scaling_settings.ForcedMosaicity(forced_mosaicity_arg);
     scaling_settings.MinCCForImage(min_image_cc / 100.0);   // --min-image-cc is in percent; the setting is a fraction
+    scaling_settings.SearchMinZeta(search_min_zeta);
     scaling_settings.OutlierRejectNsigma(
         outlier_reject_nsigma.value_or(rotation_indexing ? REJECT_OUTLIERS_DEFAULT_NSIGMA : 0.0));