diff --git a/process/JFJochProcess.cpp b/process/JFJochProcess.cpp index 27a5da27..27cc22b1 100644 --- a/process/JFJochProcess.cpp +++ b/process/JFJochProcess.cpp @@ -28,6 +28,7 @@ #include "../image_analysis/image_preprocessing/ImagePreprocessorBuffer.h" #include "../image_analysis/scale_merge/Merge.h" #include "../image_analysis/scale_merge/GlobalScale.h" +#include "../image_analysis/scale_merge/ScaleOnTheFly.h" #include "../image_analysis/scale_merge/SearchSpaceGroup.h" #include "../image_analysis/scale_merge/Combine3D.h" #include "../image_analysis/WriteReflections.h" @@ -118,6 +119,27 @@ namespace { } logger.Info("Global scaling complete ({} images)", outcomes.size()); } + + // XDS-order scaling. The rot3d combine emits fulls with partiality == 1 (image_scale_corr == 1), + // so they were only ever scaled as per-frame *partials* upstream - their per-frame scale is + // entangled with the rocking-curve/partiality model. This refits a per-frame scale directly on + // the complete reflections with the Unity model (no partiality term, G*Itrue - I_full), the way + // XDS/DIALS scale 3D-integrated fulls. A pure post-correction: it updates image_scale_corr on the + // fulls (1 -> 1/G) without re-combining. + void ScaleFulls(const DiffractionExperiment &experiment, + std::vector &fulls, int scaling_iter, + size_t nthreads, Logger &logger) { + DiffractionExperiment unity = experiment; + ScalingSettings ss = unity.GetScalingSettings(); + ss.SetPartialityModel(PartialityModel::Unity); + unity.ImportScalingSettings(ss); + + for (int i = 0; i < scaling_iter; i++) { + const auto reference = MergeAll(unity, fulls, false); + ScaleOnTheFly(unity, reference).Scale(fulls, nthreads); + } + logger.Info("Scaled fulls (XDS order, Unity model)"); + } } JFJochProcess::JFJochProcess(JFJochHDF5Reader &reader, DiffractionExperiment experiment, @@ -464,6 +486,8 @@ ProcessResult JFJochProcess::Run(JFJochProcessObserver *observer) { if (rot3d) combined = CombineRotationObservations(indexer->GetIntegrationOutcome(), experiment_, &logger, config_.observation_dump_path); + if (rot3d && config_.scale_fulls) + ScaleFulls(experiment_, combined, static_cast(config_.scaling_iter), config_.nthreads, logger); const std::vector &merge_input = rot3d ? combined : indexer->GetIntegrationOutcome(); diff --git a/process/JFJochProcess.h b/process/JFJochProcess.h index d0634421..4dafbb4a 100644 --- a/process/JFJochProcess.h +++ b/process/JFJochProcess.h @@ -56,6 +56,10 @@ struct ProcessConfig { // Global (joint) scaling instead of the alternating merge/scale loop: refine all per-image // scales, mosaicities and the shared per-HKL intensities together in one Ceres problem. bool global_scaling = false; + // XDS-order scaling: after the rot3d combine, refit a per-frame scale on the combined fulls + // (Unity model, no partiality term) so the per-frame scale is decoupled from rocking-curve + // model error. Only meaningful with -P rot3d. + bool scale_fulls = false; std::vector reference_data; // Diagnostic: if set, the -P rot3d combine writes the unmerged fulls here (for comparison vs XDS). diff --git a/tools/jfjoch_process.cpp b/tools/jfjoch_process.cpp index d387c862..54e9a025 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 << " --global-scale Joint global scaling (all images + Itrue in one fit) instead of the alternating loop; implies -M" << 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 << " -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; @@ -93,7 +94,8 @@ enum { OPT_REFERENCE_COLUMN, OPT_DUMP_OBSERVATIONS, OPT_INTEGRATOR, - OPT_GLOBAL_SCALE + OPT_GLOBAL_SCALE, + OPT_SCALE_FULLS }; static option long_options[] = { @@ -115,6 +117,7 @@ static option long_options[] = { {"wedge", optional_argument, nullptr, 'w'}, {"scale-merge", no_argument, nullptr, 'M'}, {"global-scale", no_argument, nullptr, OPT_GLOBAL_SCALE}, + {"scale-fulls", no_argument, nullptr, OPT_SCALE_FULLS}, {"refine", required_argument, nullptr, 'r'}, {"two-pass-rotation", optional_argument, nullptr, 'R'}, @@ -275,6 +278,7 @@ int main(int argc, char **argv) { std::optional rotation_indexing_range; bool run_scaling = false; bool global_scaling = false; + bool scale_fulls = false; bool anomalous_mode = false; std::optional space_group_number; std::optional fixed_reference_unit_cell; @@ -508,6 +512,10 @@ int main(int argc, char **argv) { run_scaling = true; global_scaling = true; break; + case OPT_SCALE_FULLS: + run_scaling = true; + scale_fulls = true; + break; case OPT_MIN_PARTIALITY: min_partiality = std::stod(optarg); break; @@ -791,6 +799,7 @@ int main(int argc, char **argv) { config.forced_rotation_lattice = forced_rotation_lattice; config.run_scaling = run_scaling; config.global_scaling = global_scaling; + config.scale_fulls = scale_fulls; config.scaling_iter = scaling_iter; config.reference_data = reference_data; config.observation_dump_path = dump_observations;