diff --git a/common/ScalingSettings.cpp b/common/ScalingSettings.cpp index dfda54d7..51dca65a 100644 --- a/common/ScalingSettings.cpp +++ b/common/ScalingSettings.cpp @@ -149,10 +149,34 @@ bool ScalingSettings::GetScaleFulls() const { return scale_fulls; } +ScalingSettings &ScalingSettings::SmoothGWindow(int input) { + if (input < 0) + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "Smooth-G window must be non-negative"); + smooth_g_window = input; + return *this; +} + +int ScalingSettings::GetSmoothGWindow() const { + return smooth_g_window; +} + double ScalingSettings::GetMinPartiality() const { return min_partiality; } +ScalingSettings &ScalingSettings::ForcedMosaicity(std::optional input) { + if (input.has_value() && (input.value() < GetMinMosaicity() || input.value() > GetMaxMosaicity())) + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, + "Forced mosaicity must be between " + std::to_string(GetMinMosaicity()) + + " and " + std::to_string(GetMaxMosaicity())); + forced_mosaicity = input; + return *this; +} + +std::optional ScalingSettings::GetForcedMosaicity() const { + return forced_mosaicity; +} + ScalingSettings &ScalingSettings::CaptureUncertaintyCoeff(double input) { if (input < 0.0) throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, diff --git a/common/ScalingSettings.h b/common/ScalingSettings.h index e8b50ff5..19871f3a 100644 --- a/common/ScalingSettings.h +++ b/common/ScalingSettings.h @@ -21,6 +21,7 @@ class ScalingSettings { bool merge_friedel = true; std::optional high_resolution_limit_A; std::optional wedge_for_scaling; + std::optional forced_mosaicity; // diagnostic: fix the scaling mosaicity (deg) instead of the per-image seed double min_partiality = 0.02; // Capture-aware systematic uncertainty for the rot3d combine: a full reconstructed from only // a fraction f<1 of its rocking curve is extrapolated, and the unobserved (1-f) carries a @@ -38,6 +39,11 @@ class ScalingSettings { // order). A no-op without rot3d. bool scale_fulls = false; + // Smooth the per-frame scale G across frames (centered moving average of log G over this odd + // window) before the rot3d combine, so a rocking event's partials share a consistent scale. + // 0 = off. A no-op without rot3d. + int smooth_g_window = 0; + double rfree_fraction = 0.05; IntensityFormat intensity_format = IntensityFormat::MTZ; @@ -51,11 +57,13 @@ public: ScalingSettings& HighResolutionLimit_A(double limit); ScalingSettings& HighResolutionLimit_A(std::optional limit); // nullopt clears the limit ScalingSettings& MinPartiality(double min_partiality); + ScalingSettings& ForcedMosaicity(std::optional input); ScalingSettings& CaptureUncertaintyCoeff(double input); ScalingSettings& MinCCForImage(double min_cc_for_image); ScalingSettings& OutlierRejectNsigma(double input); ScalingSettings& Combine3D(bool input); ScalingSettings& ScaleFulls(bool input); + ScalingSettings& SmoothGWindow(int input); ScalingSettings& RfreeFraction(double input); ScalingSettings& FileFormat(IntensityFormat input); @@ -81,11 +89,13 @@ public: [[nodiscard]] std::optional GetHighResolutionLimit_A() const; [[nodiscard]] double GetMinPartiality() const; + [[nodiscard]] std::optional GetForcedMosaicity() const; [[nodiscard]] double GetCaptureUncertaintyCoeff() const; [[nodiscard]] double GetMinCCForImage() const; [[nodiscard]] double GetOutlierRejectNsigma() const; [[nodiscard]] bool GetCombine3D() const; [[nodiscard]] bool GetScaleFulls() const; + [[nodiscard]] int GetSmoothGWindow() const; [[nodiscard]] double GetRfreeFraction() const; [[nodiscard]] IntensityFormat GetFileFormat() const; diff --git a/image_analysis/scale_merge/ScaleOnTheFly.cpp b/image_analysis/scale_merge/ScaleOnTheFly.cpp index 98b08a7a..3d431bd2 100644 --- a/image_analysis/scale_merge/ScaleOnTheFly.cpp +++ b/image_analysis/scale_merge/ScaleOnTheFly.cpp @@ -221,6 +221,8 @@ void ScaleOnTheFly::Scale(IntegrationOutcome &integration_outcome) const { result.mos = *integration_outcome.mosaicity_deg; else result.mos = s.GetDefaultMosaicity(); + if (const auto forced = s.GetForcedMosaicity(); forced.has_value()) + result.mos = *forced; result.wedge = rot_wedge_deg.value_or(0.0); } else { result.mos = NAN; diff --git a/process/JFJochProcess.cpp b/process/JFJochProcess.cpp index 8b2c46c9..3b5ab6fb 100644 --- a/process/JFJochProcess.cpp +++ b/process/JFJochProcess.cpp @@ -78,6 +78,44 @@ namespace { } logger.Info("Scaled fulls (XDS order, Unity model)"); } + + // Smooth the per-frame scale G across frames before the rot3d combine. ScaleOnTheFly fits each + // frame's G independently, so the few partials of one rocking event get inconsistent scales when + // weight-summed; a centered moving average of log(G) over a small odd frame window removes that + // jitter. Only G changed, so each reflection's image_scale_corr (proportional to 1/G) is rescaled + // by g_old/g_smooth. Frames without a fitted G (n &outcomes, int window, Logger &logger) { + const int n = static_cast(outcomes.size()); + const int half = window / 2; + std::vector g_smooth(n, NAN); + for (int o = 0; o < n; o++) { + double sum_log = 0.0; + int count = 0; + for (int j = std::max(0, o - half); j <= std::min(n - 1, o + half); j++) { + const auto &g = outcomes[j].image_scale_g; + if (g && std::isfinite(*g) && *g > 0.0) { + sum_log += std::log(*g); + count++; + } + } + if (count > 0) + g_smooth[o] = std::exp(sum_log / count); + } + size_t n_smoothed = 0; + for (int o = 0; o < n; o++) { + const auto &g_old = outcomes[o].image_scale_g; + if (!g_old || !std::isfinite(*g_old) || *g_old <= 0.0 || !std::isfinite(g_smooth[o])) + continue; + const double factor = *g_old / g_smooth[o]; + for (auto &r : outcomes[o].reflections) + if (std::isfinite(r.image_scale_corr)) + r.image_scale_corr = static_cast(r.image_scale_corr * factor); + outcomes[o].image_scale_g = g_smooth[o]; + n_smoothed++; + } + logger.Info("Smoothed per-frame scale G over a {}-frame window ({} of {} frames)", + window, n_smoothed, n); + } } JFJochProcess::JFJochProcess(JFJochHDF5Reader &reader, DiffractionExperiment experiment, @@ -436,6 +474,11 @@ ProcessResult JFJochProcess::Run(JFJochProcessObserver *observer) { // merging, so the error model sees counting statistics (high ISa) instead of // rocking-curve slicing scatter. const bool rot3d = experiment_.GetScalingSettings().GetCombine3D(); + if (rot3d && experiment_.GetScalingSettings().GetSmoothGWindow() > 0) { + phase("Smoothing per-frame scale G"); + SmoothImageScaleG(indexer->GetIntegrationOutcome(), + experiment_.GetScalingSettings().GetSmoothGWindow(), logger); + } std::vector combined; if (rot3d) { phase("Combining 3D partials"); diff --git a/tools/jfjoch_process.cpp b/tools/jfjoch_process.cpp index 38c4085e..9cee0d60 100644 --- a/tools/jfjoch_process.cpp +++ b/tools/jfjoch_process.cpp @@ -50,6 +50,7 @@ void print_usage() { std::cout << " Scaling and merging" << std::endl; std::cout << " -M, --scale-merge Scale and merge (refine mosaicity) and write scaled.hkl + image.dat" << std::endl; std::cout << " --scale-fulls After -P rot3d combine, refit a per-frame scale on the fulls (XDS order, Unity model); implies -M" << std::endl; + std::cout << " --smooth-g[=num] rot3d: smooth per-frame scale G over a num-frame window before the combine (default: 9 for rot3d; 0 = off)" << std::endl; std::cout << " -P, --partiality Partiality model fixed|rot|rot3d|unity (default: fixed). rot3d = rot + 3D combine of per-frame partials" << std::endl; std::cout << " -A, --anomalous Anomalous mode (don't merge Friedel pairs)" << std::endl; std::cout << " -B, --refine-bfactor Refine per image B-factor" << std::endl; @@ -57,6 +58,7 @@ void print_usage() { std::cout << " --scaling-high-resolution High resolution limit for spot finding (default: no limit)" << std::endl; std::cout << " --min-partiality Minimum partiality to accept reflection (default: 0.02)" << std::endl; std::cout << " --capture-uncertainty rot3d: systematic sigma ~num*(1-captured_fraction)*I on under-captured fulls (default: 1.0 for rot3d, 0 otherwise)" << 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: off; e.g. 6, XDS/DIALS-style)" << std::endl; std::cout << " --reject-delta-cchalf Per-crystal CC1/2-delta rejection: drop images with deltaCChalf below mean - N*stddev (default: off; e.g. 2.5)" << std::endl; std::cout << " --min-image-cc Per-image CC limit in percent (default: no limit)" << std::endl; @@ -93,7 +95,9 @@ enum { OPT_DUMP_OBSERVATIONS, OPT_INTEGRATOR, OPT_SCALE_FULLS, - OPT_CAPTURE_UNCERTAINTY + OPT_CAPTURE_UNCERTAINTY, + OPT_MOSAICITY, + OPT_SMOOTH_G }; static option long_options[] = { @@ -115,6 +119,7 @@ static option long_options[] = { {"wedge", optional_argument, nullptr, 'w'}, {"scale-merge", no_argument, nullptr, 'M'}, {"scale-fulls", no_argument, nullptr, OPT_SCALE_FULLS}, + {"smooth-g", optional_argument, nullptr, OPT_SMOOTH_G}, {"refine", required_argument, nullptr, 'r'}, {"two-pass-rotation", optional_argument, nullptr, 'R'}, @@ -129,6 +134,7 @@ static option long_options[] = { {"max-spots", required_argument, nullptr, OPT_MAX_SPOTS}, {"min-partiality", required_argument, nullptr, OPT_MIN_PARTIALITY}, {"capture-uncertainty", required_argument, nullptr, OPT_CAPTURE_UNCERTAINTY}, + {"mosaicity", required_argument, nullptr, OPT_MOSAICITY}, {"min-image-cc", required_argument, nullptr, OPT_MIN_IMAGE_CC}, {"scaling-iterations", required_argument, nullptr, OPT_SCALING_ITERATIONS}, {"scaling-high-resolution", required_argument, nullptr, OPT_SCALING_HIGH_RESOLUTION}, @@ -275,6 +281,7 @@ int main(int argc, char **argv) { std::optional rotation_indexing_range; bool run_scaling = false; bool scale_fulls = false; + std::optional smooth_g_window_arg; // --smooth-g[=window]; default 9 for rot3d, 0 (off) otherwise bool anomalous_mode = false; std::optional space_group_number; std::optional fixed_reference_unit_cell; @@ -289,6 +296,7 @@ int main(int argc, char **argv) { std::string dump_observations; // diagnostic: dump unmerged -P rot3d fulls to this path double min_partiality = 0.02; 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; int64_t scaling_iter = 3; std::optional forced_rotation_lattice; @@ -506,12 +514,18 @@ int main(int argc, char **argv) { run_scaling = true; scale_fulls = true; break; + case OPT_SMOOTH_G: + smooth_g_window_arg = optarg ? std::stoi(optarg) : 9; + break; case OPT_MIN_PARTIALITY: min_partiality = std::stod(optarg); break; case OPT_CAPTURE_UNCERTAINTY: capture_uncertainty_arg = std::stod(optarg); break; + case OPT_MOSAICITY: + forced_mosaicity_arg = std::stod(optarg); + break; case OPT_INTEGRATION_RADIUS: integration_radius_arg = optarg; break; @@ -710,6 +724,7 @@ int main(int argc, char **argv) { scaling_settings.SetPartialityModel(partiality_model); scaling_settings.Combine3D(combine_3d); scaling_settings.ScaleFulls(scale_fulls); + scaling_settings.SmoothGWindow(smooth_g_window_arg.value_or(combine_3d ? 9 : 0)); if (d_min_scale_merge) scaling_settings.HighResolutionLimit_A(d_min_scale_merge.value()); scaling_settings.MergeFriedel(!anomalous_mode); @@ -722,6 +737,7 @@ int main(int argc, char **argv) { // over-extrapolated under-captured fulls and, with the mosaicity fix, lifts rotation ISa/anomalous // substantially. Off for non-rot3d (no combine). An explicit --capture-uncertainty always wins. scaling_settings.CaptureUncertaintyCoeff(capture_uncertainty_arg.value_or(combine_3d ? 1.0 : 0.0)); + scaling_settings.ForcedMosaicity(forced_mosaicity_arg); scaling_settings.MinCCForImage(min_image_cc); if (outlier_reject_nsigma) scaling_settings.OutlierRejectNsigma(*outlier_reject_nsigma);