diff --git a/image_analysis/WriteReflections.cpp b/image_analysis/WriteReflections.cpp index dff8b25a5..bd654a999 100644 --- a/image_analysis/WriteReflections.cpp +++ b/image_analysis/WriteReflections.cpp @@ -736,12 +736,11 @@ std::vector SumRockingEvents(const std::vector & } // namespace -void WriteUnmergedMtzReflections(const std::vector &outcomes, - const UnitCell &unitCell, - const DiffractionExperiment &experiment, - bool sum_partials, - const std::string &filename, - size_t nthreads) { +gemmi::Mtz UnmergedMtz(const std::vector &outcomes, + const UnitCell &unitCell, + const DiffractionExperiment &experiment, + bool sum_partials, + size_t nthreads) { gemmi::Mtz mtz; mtz.spacegroup = &experiment.GetSpaceGroupOrP1(); mtz.set_cell_for_all(unitCell); @@ -885,13 +884,11 @@ void WriteUnmergedMtzReflections(const std::vector &outcomes // stood on that image. Turn the first indexed one back by its own angle to get the orientation of // the sweep, which is the one matrix POINTLESS also writes into every batch. std::optional lattice_at_zero; - std::optional mosaicity_deg; for (const auto &outcome : outcomes) { if (outcome.reflections.empty() || outcome.latt.CalcVolume() <= 1.0f) continue; const float mid_deg = phi_start_deg(outcome.reflections.front().image_number) + wedge_deg / 2.0f; lattice_at_zero = gon ? outcome.latt.Multiply(gon->GetTransformationAngle(mid_deg)) : outcome.latt; - mosaicity_deg = outcome.mosaicity_deg; break; } @@ -930,7 +927,6 @@ void WriteUnmergedMtzReflections(const std::vector &outcomes batch.floats[8 + 3 * i] = u[i] * z_cam; } } - batch.floats[21] = mosaicity_deg.value_or(0.0f); // crydat(0), the reflecting range batch.floats[40] = 1.0f; // scanax = [0, 0, 1]: the rotation axis IS z in the Cambridge frame batch.floats[47] = wedge_deg; batch.floats[61] = 1.0f; // e1 = scanax, the only goniostat axis @@ -975,6 +971,30 @@ void WriteUnmergedMtzReflections(const std::vector &outcomes mtz.data.swap(sorted); mtz.sort_order = {{1, 2, 3, 4, 5}}; } + return mtz; +} + +void SetUnmergedMtzMosaicity(gemmi::Mtz &mtz, const std::vector &outcomes) { + // The reflecting range of the same image the orientation matrix is taken from. + std::optional mosaicity_deg; + for (const auto &outcome : outcomes) { + if (outcome.reflections.empty() || outcome.latt.CalcVolume() <= 1.0f) + continue; + mosaicity_deg = outcome.mosaicity_deg; + break; + } + for (auto &batch : mtz.batches) + batch.floats[21] = mosaicity_deg.value_or(0.0f); // crydat(0), the reflecting range +} + +void WriteUnmergedMtzReflections(const std::vector &outcomes, + const UnitCell &unitCell, + const DiffractionExperiment &experiment, + bool sum_partials, + const std::string &filename, + size_t nthreads) { + gemmi::Mtz mtz = UnmergedMtz(outcomes, unitCell, experiment, sum_partials, nthreads); + SetUnmergedMtzMosaicity(mtz, outcomes); mtz.write_to_file(filename); } diff --git a/image_analysis/WriteReflections.h b/image_analysis/WriteReflections.h index 9aa30e89c..b44722a72 100644 --- a/image_analysis/WriteReflections.h +++ b/image_analysis/WriteReflections.h @@ -11,6 +11,8 @@ #include "../common/DiffractionExperiment.h" #include "IntegrationOutcome.h" +#include + struct MergeStatistics; struct TwinningAnalysisResult; @@ -68,6 +70,16 @@ void WriteUnmergedMtzReflections(const std::vector &outcomes const std::string &filename, size_t nthreads = 1); +// The same file in two steps, for a caller that builds it while the outcomes' per-image scale fields +// are still being written: UnmergedMtz reads everything but mosaicity_deg, which +// SetUnmergedMtzMosaicity then copies into the batch headers. +gemmi::Mtz UnmergedMtz(const std::vector &outcomes, + const UnitCell &unitCell, + const DiffractionExperiment &experiment, + bool sum_partials, + size_t nthreads = 1); +void SetUnmergedMtzMosaicity(gemmi::Mtz &mtz, const std::vector &outcomes); + void WriteReflections(const std::vector &reflections, const UnitCell &unitCell, const DiffractionExperiment &experiment, diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index 460c308fd..ba9b26063 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -6104,6 +6104,9 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b end_msg.unit_cell = result.consensus_cell; } + // The unmerged MTZ, when it is built beside the P1 cross-check merge (see there). + std::future unmerged_mtz; + // Scaling and merging (full analysis only). if (full && !cancelled_ && (result.indexing_rate.value_or(0.0f) > 0 || result.rotation_lattice_found) && (config_.run_scaling || !config_.reference_data.empty())) { @@ -8742,6 +8745,17 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b const double em_a = result.error_model_a; const double em_b = result.error_model_b; const auto res_fit = result.resolution_fit_A; + // The unmerged MTZ below does not depend on this merge, so it is built meanwhile, from + // the experiment as it stands in the determined group. The merge rewrites each image's + // mosaicity, which the file's batch headers carry, so that is filled in only after it. + if (full && write_files && !geometry_prepass && !superseded && result.consensus_cell + && config_.export_unmerged) + unmerged_mtz = std::async(std::launch::async, + [&outcomes = indexer->GetIntegrationOutcome(), + cell = *result.consensus_cell, x = experiment_, + nthreads = config_.nthreads] { + return UnmergedMtz(outcomes, cell, x, true, nthreads); + }); // Both the merge and the MTZ read the group from the experiment, so it is set for // the whole of it and restored after. experiment_.SpaceGroupNumber(1); @@ -8841,8 +8855,13 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b if (config_.export_unmerged) { if (observer) observer->OnPhase("Writing unmerged reflections"); const std::string path = config_.output_prefix + "_unmerged.mtz"; - WriteUnmergedMtzReflections(indexer->GetIntegrationOutcome(), *result.consensus_cell, - experiment_, true, path, config_.nthreads); + if (unmerged_mtz.valid()) { + gemmi::Mtz mtz = unmerged_mtz.get(); + SetUnmergedMtzMosaicity(mtz, indexer->GetIntegrationOutcome()); + mtz.write_to_file(path); + } else + WriteUnmergedMtzReflections(indexer->GetIntegrationOutcome(), *result.consensus_cell, + experiment_, true, path, config_.nthreads); logger.Info("Unmerged observations written to {}", path); } if (config_.export_unmerged_partials) {