Enable reading MTZ files with reference Fcalc

This commit is contained in:
2026-05-13 14:22:39 +02:00
parent 08bf186766
commit 6aa6890a8b
14 changed files with 100 additions and 20 deletions
+23 -5
View File
@@ -25,6 +25,7 @@
#include "../image_analysis/IndexAndRefine.h"
#include "../receiver/JFJochReceiverPlots.h"
#include "../compression/JFJochCompressor.h"
#include "../image_analysis/LoadFCalcFromMtz.h"
#include "../image_analysis/scale_merge/FrenchWilson.h"
#include "../image_analysis/scale_merge/SearchSpaceGroup.h"
#include "../image_analysis/WriteMmcif.h"
@@ -149,6 +150,7 @@ int main(int argc, char **argv) {
bool refine_bfactor = false;
bool refine_wedge = false;
std::optional<double> wedge_for_scaling;
std::string ref_mtz;
IndexingAlgorithmEnum indexing_algorithm = IndexingAlgorithmEnum::Auto;
@@ -163,11 +165,14 @@ int main(int argc, char **argv) {
}
int opt;
while ((opt = getopt(argc, argv, "o:N:s:e:vc:R::FX:xd:S:MP:AD:C:T:t:Bw::")) != -1) {
while ((opt = getopt(argc, argv, "o:N:s:e:vc:R::FX:xd:S:MP:AD:C:T:t:Bw::z:")) != -1) {
switch (opt) {
case 'o':
output_prefix = optarg;
break;
case 'z':
ref_mtz = optarg;
break;
case 'N':
nthreads = atoi(optarg);
break;
@@ -319,6 +324,12 @@ int main(int argc, char **argv) {
logger.Info("Loaded dataset from {}", input_file);
std::vector<MergedReflection> reference_data;
if (!ref_mtz.empty()) {
reference_data = LoadFCalcFromMtz(ref_mtz);
logger.Info("Loaded {} reflections from {} MTZ file", reference_data.size(), ref_mtz);
}
// 2. Setup Experiment & Components
DiffractionExperiment experiment(dataset->experiment);
experiment.BitDepthImage(32).Compression(CompressionAlgorithm::BSHUF_LZ4);
@@ -579,10 +590,17 @@ int main(int argc, char **argv) {
auto scale_start = std::chrono::steady_clock::now();
for (int i = 0; i < 3; i++) {
auto iter_start = std::chrono::steady_clock::now();
auto scale_result = indexer.ScaleAllImages();
end_msg.image_scale_factor = scale_result.image_scale_g;
scale_result.SaveToFile(output_prefix + "_iter" + std::to_string(i) + "_scale.dat");
if (reference_data.empty()) {
auto merge_result = indexer.Merge(false);
auto scale_result = indexer.ScaleAllImages(merge_result.merged);
end_msg.image_scale_factor = scale_result.image_scale_g;
scale_result.SaveToFile(output_prefix + "_iter" + std::to_string(i) + "_scale.dat");
} else {
auto scale_result = indexer.ScaleAllImages(reference_data);
end_msg.image_scale_factor = scale_result.image_scale_g;
scale_result.SaveToFile(output_prefix + "_iter" + std::to_string(i) + "_scale.dat");
}
auto iter_end = std::chrono::steady_clock::now();
double iter_time = std::chrono::duration<double>(iter_end - iter_start).count();
@@ -592,7 +610,7 @@ int main(int argc, char **argv) {
double scale_time = std::chrono::duration<double>(scale_end - scale_start).count();
auto merge_start = std::chrono::steady_clock::now();
auto merge_result = indexer.Merge();
auto merge_result = indexer.Merge(true);
auto merge_end = std::chrono::steady_clock::now();
double merge_time = std::chrono::duration<double>(merge_end - merge_start).count();