diff --git a/common/ScalingSettings.cpp b/common/ScalingSettings.cpp index 9a43f393..617ff2ac 100644 --- a/common/ScalingSettings.cpp +++ b/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"); diff --git a/common/ScalingSettings.h b/common/ScalingSettings.h index abef9dfa..60b98b0d 100644 --- a/common/ScalingSettings.h +++ b/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; diff --git a/image_analysis/scale_merge/RotationScaleMerge.cpp b/image_analysis/scale_merge/RotationScaleMerge.cpp index c8e3f4ae..9b698dbe 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.cpp +++ b/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 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 reject(n_frames, 0); diff --git a/image_analysis/scale_merge/RotationScaleMerge.h b/image_analysis/scale_merge/RotationScaleMerge.h index adb37fa9..f93ec27c 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.h +++ b/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; diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index cde3d25e..30388b68 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -109,6 +109,7 @@ void print_usage() { std::cout << " --min-captured-fraction 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 Diagnostic: fix the scaling mosaicity (deg) instead of the per-image seed" << std::endl; std::cout << " --reject-outliers 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 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 Per-image CC limit in percent (default: no limit)" << std::endl; std::cout << " --scaling-iterations Number of scaling iterations with no reference data (default: 3)" << std::endl; std::cout << " -z, --reference-mtz 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 capture_uncertainty_arg; // explicit --capture-uncertainty; default depends on rot3d std::optional 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 forced_rotation_lattice; std::optional 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));