diff --git a/tools/jfjoch_process.cpp b/tools/jfjoch_process.cpp index 591f7840..218b6a17 100644 --- a/tools/jfjoch_process.cpp +++ b/tools/jfjoch_process.cpp @@ -32,7 +32,7 @@ #include "../image_analysis/WriteReflections.h" #include "../image_analysis/UpdateReflectionResolution.h" -void print_usage(Logger &logger) { +void print_usage() { std::cout << "Usage ./jfjoch_analysis {} " << std::endl; std::cout << "Options:" << std::endl; std::cout << " -o, --output-prefix Output file prefix (default: output)" << std::endl; @@ -46,30 +46,30 @@ void print_usage(Logger &logger) { std::cout << " Spot finding" << std::endl; std::cout << " --spot-sigma Noise sigma level for spot finding (default: 3.0)" << std::endl; std::cout << " --spot-threshold Photon count threshold for spot finding (default: 10)" << std::endl; - std::cout << " --spot-resolution High resolution limit for spot finding (default: 1.5)" << std::endl; + std::cout << " --spot-high-resolution High resolution limit for spot finding (default: 1.5)" << std::endl; std::cout << " --max-spots Max spot count (default: 250)" << std::endl; std::cout << std::endl; std::cout << " Indexing" << std::endl; std::cout << " -R, --rotation-indexing[=num] Rotation indexing (optional: min angular range deg)" << std::endl; std::cout << " -X, --indexing-algorithm Indexing algorithm (FFBIDX|FFT|FFTW|Auto|None)" << std::endl; - std::cout << " -F, --fft-indexing Use FFT indexing algorithm (shortcut for -X FFT)" << std::endl; std::cout << " -S, --space-group Space group number - used for both indexing and scaling" << std::endl; std::cout << " -C, --unit-cell Fix reference unit cell: \"a,b,c,alpha,beta,gamma\"" << std::endl; std::cout << " --no-beam-center-refine No least-square beam center refinement" << std::endl; std::cout << std::endl; 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 << " -P, --partiality Partiality refinement fixed|rot|unity (default: fixed)" << std::endl; + std::cout << " -M, --scale-merge Scale and merge (refine mosaicity) and write scaled.hkl + image.dat" << std::endl; + std::cout << " -P, --partiality Partiality refinement fixed|rot|unity (default: fixed)" << 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; std::cout << " -w, --wedge[=num] Refine image wedge during scaling with starting wedge value" << std::endl; + 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 << " --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; -} + } enum { OPT_SPOT_SIGMA = 1000, @@ -79,7 +79,8 @@ enum { OPT_NO_BEAM_CENTER_REFINE, OPT_MIN_PARTIALITY, OPT_MIN_IMAGE_CC, - OPT_SCALING_ITERATIONS + OPT_SCALING_ITERATIONS, + OPT_SCALING_HIGH_RESOLUTION }; static option long_options[] = { @@ -88,7 +89,7 @@ static option long_options[] = { {"threads", required_argument, nullptr, 'N'}, {"start-image", required_argument, nullptr, 's'}, {"end-image", required_argument, nullptr, 'e'}, -{"stride", required_argument, nullptr, 't'}, + {"stride", required_argument, nullptr, 't'}, {"rotation-indexing", optional_argument, nullptr, 'R'}, {"indexing-algorithm", required_argument, nullptr, 'X'}, {"unit-cell", required_argument, nullptr, 'C'}, @@ -102,12 +103,13 @@ static option long_options[] = { {"spot-sigma", required_argument, nullptr, OPT_SPOT_SIGMA}, {"spot-threshold", required_argument, nullptr, OPT_SPOT_THRESHOLD}, - {"spot-resolution", required_argument, nullptr, OPT_SPOT_RESOLUTION}, + {"spot-high-resolution", required_argument, nullptr, OPT_SPOT_RESOLUTION}, {"max-spots", required_argument, nullptr, OPT_MAX_SPOTS}, {"no-beam-center-refine", no_argument, nullptr, OPT_NO_BEAM_CENTER_REFINE}, {"min-partiality", required_argument, nullptr, OPT_MIN_PARTIALITY}, {"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}, {nullptr, 0, nullptr, 0} }; @@ -219,7 +221,7 @@ int main(int argc, char **argv) { std::optional d_min_scale_merge; if (argc == 1) { - print_usage(logger); + print_usage(); exit(EXIT_FAILURE); } @@ -268,7 +270,7 @@ int main(int argc, char **argv) { indexing_algorithm = IndexingAlgorithmEnum::None; else { logger.Error("Invalid indexing algorithm: {}", alg); - print_usage(logger); + print_usage(); exit(EXIT_FAILURE); } break; @@ -279,7 +281,7 @@ int main(int argc, char **argv) { logger.Error( "Invalid unit cell. Expected: \"a,b,c,alpha,beta,gamma\" (6 floats, comma-separated, no spaces). Got: {}", optarg ? optarg : ""); - print_usage(logger); + print_usage(); exit(EXIT_FAILURE); } fixed_reference_unit_cell = uc; @@ -317,7 +319,7 @@ int main(int argc, char **argv) { partiality_model = PartialityModel::Rotation; else { logger.Error("Invalid partiality mode: {}", optarg); - print_usage(logger); + print_usage(); exit(EXIT_FAILURE); } break; @@ -350,6 +352,9 @@ int main(int argc, char **argv) { case OPT_MIN_IMAGE_CC: min_image_cc = std::stod(optarg); break; + case OPT_SCALING_HIGH_RESOLUTION: + d_min_scale_merge = atof(optarg); + break; case OPT_SCALING_ITERATIONS: scaling_iter = atoi(optarg); if (scaling_iter <= 0) { @@ -359,14 +364,14 @@ int main(int argc, char **argv) { break; default: - print_usage(logger); + print_usage(); exit(EXIT_FAILURE); } } if (optind != argc - 1) { logger.Error("Input file not specified"); - print_usage(logger); + print_usage(); exit(EXIT_FAILURE); } diff --git a/tools/jfjoch_scale.cpp b/tools/jfjoch_scale.cpp index 4a8c9a87..830d072c 100644 --- a/tools/jfjoch_scale.cpp +++ b/tools/jfjoch_scale.cpp @@ -11,6 +11,7 @@ #include #include #include +#include #include "../reader/JFJochHDF5Reader.h" #include "../common/Logger.h" @@ -31,27 +32,60 @@ #include "../image_analysis/WriteReflections.h" #include "../image_analysis/UpdateReflectionResolution.h" -void print_usage(Logger &logger) { - logger.Info("Usage ./jfjoch_analysis {} "); - logger.Info("Options:"); - logger.Info(" -o Output file prefix (default: output)"); - logger.Info(" -N Number of threads (default: 1)"); - logger.Info(" -s Start image number (default: 0)"); - logger.Info(" -e End image number (default: all)"); - logger.Info(" -v Verbose output"); - logger.Info(""); - logger.Info(" Scaling and merging"); - logger.Info(" -S Space group number"); - logger.Info(" -D High resolution limit for scaling/merging (default: 0.0; no limit)"); - logger.Info(" -P Partiality refinement fixed|rot|unity (default: fixed)"); - logger.Info(" -A Anomalous mode (don't merge Friedel pairs)"); - logger.Info(" -B Refine per image B-factor"); - logger.Info(" -w Refine image wedge during scaling with starting wedge value"); - logger.Info(" -p Minimum partiality to accept reflection (default: 0.02)"); - logger.Info(" -q Per-image CC limit in percent (default: no limit)"); - logger.Info(" -i Number of scaling iterations with no reference data (default: 3) "); +void print_usage() { + std::cout << "Usage ./jfjoch_scale {} " << std::endl; + std::cout << "Options:" << std::endl; + std::cout << " -o, --output-prefix Output file prefix (default: output)" << std::endl; + std::cout << " -N, --threads Number of threads (default: 1)" << std::endl; + std::cout << " -s, --start-image Start image number (default: 0)" << std::endl; + std::cout << " -e, --end-image End image number (default: all)" << std::endl; + std::cout << " -v, --verbose Verbose output" << std::endl; + std::cout << "" << std::endl; + std::cout << " Scaling and merging" << std::endl; + std::cout << " -S, --space-group Space group number" << std::endl; + std::cout << + " -D, --scaling-high-resolution High resolution limit for scaling/merging (default: 0.0; no limit)" + << std::endl; + std::cout << " -P, --partiality Partiality refinement fixed|rot|unity (default: fixed)" << + 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; + std::cout << " -w, --wedge[=num] Refine image wedge during scaling with starting wedge value" + << std::endl; + std::cout << " --min-partiality Minimum partiality to accept reflection (default: 0.02)" << + 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; } +enum { + OPT_MIN_PARTIALITY = 1000, + OPT_MIN_IMAGE_CC, + OPT_SCALING_ITERATIONS +}; + + +static option long_options[] = { + {"verbose", no_argument, nullptr, 'v'}, + {"output-prefix", required_argument, nullptr, 'o'}, + {"threads", required_argument, nullptr, 'N'}, + {"start-image", required_argument, nullptr, 's'}, + {"end-image", required_argument, nullptr, 'e'}, + {"reference-mtz", required_argument, nullptr, 'z'}, + {"space-group", required_argument, nullptr, 'S'}, + {"partiality", required_argument, nullptr, 'P'}, + {"anomalous", no_argument, nullptr, 'A'}, + {"refine-bfactor", no_argument, nullptr, 'B'}, + {"wedge", optional_argument, nullptr, 'w'}, + {"scaling-high-resolution", required_argument, nullptr, 'D'}, + {"min-partiality", required_argument, nullptr, OPT_MIN_PARTIALITY}, + {"min-image-cc", required_argument, nullptr, OPT_MIN_IMAGE_CC}, + {"scaling-iterations", required_argument, nullptr, OPT_SCALING_ITERATIONS}, + {nullptr, 0, nullptr, 0} +}; + int main(int argc, char **argv) { RegisterHDF5Filter(); @@ -80,18 +114,21 @@ int main(int argc, char **argv) { std::optional d_min_scale_merge; if (argc == 1) { - print_usage(logger); + print_usage(); exit(EXIT_FAILURE); } int opt; - while ((opt = getopt(argc, argv, "o:z:N:s:e:p:q:Bw::vD:S:A:P:")) != -1) { + int option_index = 0; + const char *short_opts = "vo:N:s:e:z:S:P:ABw::D:p:q:i:"; + + while ((opt = getopt_long(argc, argv, short_opts, long_options, &option_index)) != -1) { switch (opt) { case 'o': output_prefix = optarg; break; - case 'z': - ref_mtz = optarg; + case 'v': + verbose = true; break; case 'N': nthreads = atoi(optarg); @@ -102,33 +139,12 @@ int main(int argc, char **argv) { case 'e': end_image = atoi(optarg); break; - case 'p': - min_partiality = std::stod(optarg); - break; - case 'q': - min_image_cc = std::stod(optarg); - break; - case 'B': - refine_bfactor = true; - break; - case 'w': - refine_wedge = true; - if (optarg) - wedge_for_scaling = std::stod(optarg); - break; - case 'v': - verbose = true; - break; - case 'D': - d_min_scale_merge = atof(optarg); - logger.Info("High resolution limit for scaling/merging set to {:.2f} A", d_min_scale_merge.value()); + case 'z': + ref_mtz = optarg; break; case 'S': space_group_number = atoi(optarg); break; - case 'A': - anomalous_mode = true; - break; case 'P': if (strcmp(optarg, "unity") == 0) partiality_model = PartialityModel::Unity; @@ -138,26 +154,44 @@ int main(int argc, char **argv) { partiality_model = PartialityModel::Rotation; else { logger.Error("Invalid partiality mode: {}", optarg); - print_usage(logger); + print_usage(); exit(EXIT_FAILURE); } break; - case 'i': + case 'A': + anomalous_mode = true; + break; + case 'B': + refine_bfactor = true; + break; + case 'w': + refine_wedge = true; + if (optarg) + wedge_for_scaling = std::stod(optarg); + break; + case 'D': + d_min_scale_merge = atof(optarg); + logger.Info("High resolution limit for scaling/merging set to {:.2f} A", d_min_scale_merge.value()); + break; + case OPT_MIN_PARTIALITY: + min_partiality = std::stod(optarg); + break; + case OPT_MIN_IMAGE_CC: + min_image_cc = std::stod(optarg); + break; + case OPT_SCALING_ITERATIONS: scaling_iter = atoi(optarg); if (scaling_iter <= 0) { logger.Error("Invalid scaling iteration count: {}", scaling_iter); exit(EXIT_FAILURE); } break; - default: - print_usage(logger); - exit(EXIT_FAILURE); } } if (optind != argc - 1) { logger.Error("Input file not specified"); - print_usage(logger); + print_usage(); exit(EXIT_FAILURE); } @@ -252,7 +286,7 @@ int main(int argc, char **argv) { auto refl_stats = UpdateReflectionResolution(experiment.GetUnitCell().value(), reflections); logger.Info("Read {} reflections from {} images", refl_stats.n_reflections, refl_stats.n_images); - + std::vector mosaicity(end_image - start_image + 1); for (int i = 0; i < end_image - start_image + 1; i++) { mosaicity[i] = dataset->mosaicity_deg[start_image + i]; @@ -270,7 +304,7 @@ int main(int argc, char **argv) { scale_result = scaling.Scale(reflections, mosaicity, nthreads); } else { ScaleOnTheFly scaling(experiment, reference_data); - scale_result = scaling.Scale(reflections,mosaicity, nthreads); + scale_result = scaling.Scale(reflections, mosaicity, nthreads); } scale_result.SaveToFile(output_prefix + "_iter" + std::to_string(i) + "_scale.dat"); @@ -292,7 +326,7 @@ int main(int argc, char **argv) { auto merge_result = MergeAll(experiment, reflections, merge_mask); auto merge_stats = MergeStats(experiment, merge_result, reflections, experiment.GetUnitCell().value(), merge_mask, - reference_data); + reference_data); auto merge_end = std::chrono::steady_clock::now(); double merge_time = std::chrono::duration(merge_end - merge_start).count();